Momentum flow in black-hole binaries: II. Numerical simulations of equal-mass, 

head-on mergers with antiparallel spins. 



Geoffrey Lovelace,^ Yanbei Chen,^ Michael Cohen, ^ Jeffrey D. Kaplan,^ Drew Keppel,^ 
Keith D. Matthews,^ David A. Nichols,^ Mark A. Scheel,^ and Ulrich Sperhake^ 

^ Center for Radiophysics and Space Research, Cornell University, Rhaca, New York, 14853 
^Theoretical Astrophysics 350-17, California Institute of Technology, Pasadena, CA 91125 

(Dated: July 5, 2009) 

Research on extracting science from binary-black-hole (BBH) simulations has often adopted a 
"scattering matrix" perspective: given the binary's initial parameters, what are the final hole's 
parameters and the emitted gravitational waveform? In contrast, we are using BBH simulations 
to explore the nonlinear dynamics of curved spacetime. Focusing on the head-on plunge, merger, 
and ringdown of a BBH with transverse, antiparallel spins, we explore numerically the momentum 
flow between the holes and the surrounding spacetime. We use the Landau-Lifshitz field-theory- 
in-flat-spacetime formulation of general relativity to define and compute the density of field energy 
and field momentum outside horizons and the energy and momentum contained within horizons, 
and we define the effective velocity of each apparent and event horizon as the ratio of its enclosed 
momentum to its enclosed mass-energy. We find surprisingly good agreement between the horizons' 
effective and coordinate velocities. During the plunge, the holes experience a frame-dragging-induced 
acceleration orthogonal to the plane of their spins and their infall ("downward"), and they reach 
downward speeds of order 1000 km/s. When the common apparent horizon forms (and when the 
event horizons merge and their merged neck expands) , the horizon swallows upward field momentum 
that resided between the holes, causing the merged hole to accelerate in the opposite ("upward") 
direction. As the merged hole and the field energy and momentum settle down, a pulsational burst 
of gravitational waves is emitted, and the merged hole has a final effective velocity of about 20 
km/s upward, which agrees with the recoil velocity obtained by measuring the linear momentum 
carried to infinity by the emitted gravitational radiation. To investigate the gauge dependence of 
our results, we compare pseudospectral and moving-puncture evolutions of physically similar initial 
data; although spectral and puncture simulations use different gauge conditions, we find remarkably 
good agreement for our results in these two cases. We also compare our simulations with the 
post-Newtonian trajectories and near-field energy-momentum. 

PACS numbers: 04.25.D-, 04.25.dg, 04.25.Nx, 04.70.-s, 97.60.Lf 



I. INTRODUCTION 
A. Motivation 

Following Pretorius's 2005 breakthrough [1 , several re- 
search groups have developed codes to solve Einstein's 
equations numerically for the inspiral, merger, and ring- 
down of colliding binary black holes (BBHs). Most simu- 
lations of BBH mergers to date have adopted the moving- 
puncture method [2 , 3 , and spectral methods [4 have 
also successfully simulated BBH mergers. 

A major goal of current research is to successfully ex- 
tract the physical content of these simulations. Typi- 
cally, efforts toward this goal adopt a "scattering matrix" 
approach. Information obtained from numerical simula- 
tions on a finite set of islands in the seven-dimensional^ 
parameter space is being extrapolated, by various re- 
search groups, to design complicated functions that give 



^ One parameter for the mass ratio and six for the individual spins; 
additional parameters might arise from eccentric orbits and the 
apparent dependence, in at least some configurations, of the re- 
coil on the initial phase of the binary. 



the final parameters of the merged hole and the emit- 
ted gravitational waveforms as functions of the binary's 
initial parameters. 

In this paper, however, we take a different perspec- 
tive: we focus our attention on the nonlinear dynamics 
of curved spacetime during the holes' merger and ring- 
down. Following Ref. [5] (paper I in this series), our goal 
is to develop physical insight into the behavior of highly 
dynamical spacetimes such as the strong-field region near 
the black-hole horizons in a merging binary. As in pa- 
per I, we focus this study on the distribution and flow of 
linear momentum in BBH spacetimes. In contrast to pa- 
per I's description of the pre-merger motion of the holes 
in the post-Newtonian approximation, in this paper we 
study the momentum flow during the plunge, merger, 
and ringdown of merging black holes in fully relativistic 
simulations. 



B. Linear momentum flow in BBHs and gauge 
dependence 

Typically, numerical simulations calculate only the 
total linear momentum of a BBH system and ignore 
the (gauge-dependent) linear momenta of the individual 
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black holes. However, linear momentum has been consid- 
ered by Krishnan, Lousto and Zlochower [6^. Inspired by 
the success of quasilocal angular momentum (see, e.g., [7 
for a review) as a tool for measuring the spin of an indi- 
vidual black hole, Krishnan and colleagues proposed an 
analogous (but gauge-dependent) formula for the quasilo- 
cal linear momentum, and they calculate this quasilocal 
linear momentum for, e.g., the highly-spinning, unequal- 
mass BBH simulations in Ref. |8]. This quasilocal linear 
momentum is also used to define an orbital angular mo- 
mentum in Ref. pi. 

In this paper, we adopt a different, complementary 
method for measuring the holes' linear momenta: for 
the first time, we apply the Landau-Lifshitz momentum- 
fiow formalism (described in paper I and summarized in 
Sec. [ll| to numerical simulations of merging black holes. 
In this formalism, a mapping between the curved space- 
time and an auxiliary fiat spacetime (AFS) is chosen, 
and general relativity is reinterpreted as a field theory 
defined on this fiat spacetime. The AFS has a set of 
translational Killing vectors which we use to define a lo- 
calized, conserved linear momentum. In particular, we 
calculate i) a momentum density, ii) the momentum en- 
closed by horizons, and iii) the momentum enclosed by 
distant coordinate spheres. In the asymptotically fiat re- 
gion around a source, there is a preferred way to choose 
the mapping between the curved spacetime and the AFS; 
consequently, in this limit item iii) is gauge- invariant. 
In general, though, the choice of mapping is arbitrary, 
and it follows that items i) and ii) are necessarily gauge- 
dependent. 

By examining the linear momentum fiow in a dy- 
namical spacetime — and living with the inevitable gauge 
dependence — we hope to develop strong intuition for the 
behavior of BBHs. As in paper I, we envision different 
numerical relativity groups choosing "preferred" gauges 
based on the coordinates of their numerical simulations. 
While there is no reason, a priori^ why simulations in 
different gauges should agree, one of our hopes from pa- 
per I is realized for the cases we consider; namely, in 
this paper, we calculate the horizon-enclosed momentum 
using spectral and moving-puncture evolutions of similar 
initial data, and we do find surprisingly good agreement 
(cf. Figs. Hand [15|, even though the simulations use 
manifestly different gauge conditions [Eqs. (14) for the 
spectral simulations and Eqs. (B19)-(B20) for the punc- 



ture simulations]. These are two of the most commonly 
used gauge conditions in numerical relativity. 

Therefore, we continue to hope that in general — for 
the gauges commonly used in numerical simulations — 
the momentum distributions for evolutions of physically 
similar initial data will turn out to be at least qualita- 
tively similar. If further investigation reveals this to be 
the case, then different research groups can simply use 
the coordinates used in the their simulations as the "pre- 
ferred coordinates" for constructing the mapping to the 
AFS. Otherwise, we would advocate (as in Sec. I C of 
paper I) that different numerical-relativity groups con- 



struct the mapping to the AFS by first agreeing on a 
choice of "preferred" coordinates (e.g., a particular har- 
monic gauge) and then transforming the results of their 
simulations to those coordinates. 



C. BBH mergers with recoil 

A particularly important application of this approach 
is an exploration of the momentum fiow in BBH merg- 
ers with recoil. The gravitational recoil or kick effect 
arising in a BBH coalescence has attracted a great deal 
of attention in recent years in the context of a vari- 
ety of astrophysical scenarios including the structure 
of galaxies [lOj [TTJ [12], the reionization history of the 
universe [13], the assembly of supermassive black holes 
[14, [161 HZl dB! and direct observational signatures 
[191 [ini [21]. For a long time, estimates of the re- 
coil magnitude were based on approximative techniques 
[22l [23l [24[ [25] ; accurate calculations in the framework of 
fully nonlinear general relativity have only become pos- 
sible in the aftermath of important breakthroughs in the 
field of numerical relativity [11 [2] [3] . 

Several groups have used numerical simulations to 
study the kick resulting from the merger of non-spinning 
and spinning binaries (see, e.g., [26 } [27 l [28 } [29 } [50 1 [3T] ) . 
Most remarkably, recoil velocities of several thousand 
km/s have been found for binaries with equal and op- 
posite spins in the orbital plane [30l |32l |33] , and variants 
thereof with hyperbolic orbits even generate 10^ km/s 
[34 . Given the enormous astrophysical repercussions of 
such large recoil velocities, the community is now us- 
ing various approaches to obtain a better understand- 
ing of the kick as a function of the initial BBH parame- 
ters [35] [36l [371 EHl [SH [40] resulting in phenomenological 
fitting formulas; see [H [9] [38l [411 [HI [Sj and references 
therein. 

On the other hand, our understanding of the local 
dynamics in these extraordinarily violent events is still 
rather limited. Some insight into the origin of the holes' 
kick velocity has been obtained by examining the indi- 
vidual multipole moments of the emitted gravitational 
waves UK [45] and by approximating the recoil ana- 
lytically using post-Newtonian [24l [46], effective-one- 
body [25], and black-hole-perturbation theory [47] • An 
intuitive picture describing aspects of the so-called super- 
kick configurations generating velocities in the thousands 
of km/s has been given in terms of the frame-dragging 
effect (cf. Fig. 5 of Ref. [48 ). 

Investigating the momentum distribution and fiow in 
recoiling BBH mergers could help to build further intu- 
ition into the nonlinear dynamics of the spacetime and 
their infiuence on the formation of kicks. Paper I made 
some headway into the former issue but could not address 
the latter. Specifically, paper I examined the distribution 
and the fiow of linear momentum in BBH spacetimes us- 
ing the Landau-Lifshitz formalism in the post-Newtonian 
approximation. It then specialized this approach to the 
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FIG. 1: Initial configuration of tfie fiead-on BBH considered 
in tfiis paper. Tfie fioles move primarily along the x axis, but 
they also accelerate in the —y (downward) direction due to 
frame dragging. See Table [l] for the value of d = 2xo. 



extreme-kick configuration [30l |32J [33] , which is a system 
of inspiraling, BBHs with equal and anti-parallel spins in 
the orbital plane. During inspiral, the two black holes si- 
multaneously and sinusoidally bob perpendicularly to the 
orbital plane; in paper I, this motion was first recognized 
as arising from the combined effect of frame dragging 
and spin-curvature coupling and then was found to arise 
from the exchange of momentum between the near-zone 
gravitational field and the black holes. 

Because paper I analyzed the system at a post- 
Newtonian level, its analysis could not be extended to 
merger and beyond. Consequently, it was not possible 
to address how the nonlinear dynamics in the pre-merger 
near zone transitions into the final behavior of the merged 
black hole. This paper (paper II) lets us begin to address 
this transition as we study momentum flow during the 
plunge, merger, and ringdown of BBHs in full numerical 
relativity. Our study allows us, for example, to examine 
how accurately Pretorius's intuitive picture applies dur- 
ing the merger and ringdown of a recoiling BBH merger. 



D. Overview and summary 

As a first step toward analyzing the momentum flow 
in superkicks, in this paper we apply the Landau-Lifshitz 
momentum- flow formalism to a much simpler case: the 
head-on plunge, merger, and ringdown of an equal-mass 
BBH. The holes initially have antiparallel spins of equal 
magnitude that are transverse to the holes' head-on mo- 
tion (Fig. [T]). Primarily, the holes simply fall toward 
each other in the ±x direction. However, each hole's 
spin drags the space around itself, causing the other hole 
to accelerate in the downward, —y direction. 

How does this frame dragging relate to the final kick 
velocity of the merged hole? To address this question, we 
compute the 4-momentum inside each apparent hori- 
zon using the Landau-Lifshitz formalism; we then define 
an effective velocity as 
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FIG. 2: The effective velocity v^-^ for the individual (red dot- 
ted line) and common (green dashed line) apparent horizons 
and for the event horizon (black solid line). The inset shows 
the velocity of the common apparent horizon at late times. 



In Sec. |IV[ we find that this effective velocity behaves 
similarly to the horizons' coordinate velocities. 

The effective y velocity for the spectral simulation de- 

Fig. 



scribed in Sec. HI A 2 
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IS shown in Fig. |2] Before the 
merger, the individual apparent horizons do indeed accel- 
erate in the —y ("down") direction, eventually reaching 
velocities of order 10^ km/s. However, when the common 
apparent horizon forms, it pulsates; during the first half- 
pulsation, the horizon expands and accelerates to ~ 10^ 
km/s in the up {-\-y) direction. This happens because as 
the common horizon forms and expands, it swallows not 
only the downward linear momentum inside each individ- 
ual horizon but also a large amount of upward momentum 
in the gravitational field between the holes (Fig.[3|. Dur- 
ing the next half-pulsation, as the horizon shape changes 
from oblate to prolate (cf. Fig. pT|, the horizon swallows 
a net downward momentum, thereby losing most of its 
upward velocity. Eventually, after strong damping of the 
pulsations, the common horizon settles down to a very 
small velocity of about 23 km/s in the -\-y direction (inset 
of Fig.[2|, which (Sec. [lV|) is consistent with the kick ve- 
locity inferred from the emitted gravitational radiation. 

This momentum flow between field and holes is also 
described quite beautifully in the language of the holes' 
event horizon. Unlike apparent horizons, the event hori- 
zon evolves and expands continuously in time, rather 
than discontinuously. As the event horizon expands, it 
continuously swallows surrounding field momentum, and 
that swallowing produces a continuous evolution of the 
event horizon's velocity, an evolution that is nearly the 
same as for the apparent-horizon velocity. Figure |2] shows 
how the effective velocity of the event horizon smoothly 
transitions from matching the individual apparent hori- 
zons' velocities to matching the common apparent hori- 
zon's velocity. For further details, see Sec. |IV A2| and 
especially Figs. [13] and [14] 

In the remainder of this paper, we discuss our results 
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Time = 34.73 Madm 




FIG. 3: A contour plot of the y component of the momentum 
density at the moment when the common apparent horizon 
forms. The common horizon encloses the momentum inside 
the individual horizons and also the momentum in the gravi- 
tational field. The grey-shaded region and solid, red contours 
indicate positive momentum density, while the white-shaded 
region and blue, dashed contours indicate negative momen- 
tum density. The individual apparent horizons are shaded 
black, and the common apparent horizon is shown as a thick 
black line. 



and the simulations that are used to obtain them. In 
Sec. [nj we briefly review the Landau-Lifshitz formalism 
and momentum conservation. The simulations them- 
selves are presented in Sec. |III[ W e analyze the simula- 
tions' momentum flow in Sec. lIVI and conclude in Sec. IVl 
In the appendices, we describe in greater depth the nu- 
merical methods used for the simulations presented in 
this paper. 



II. 4-MOMENTUM CONSERVATION IN THE 
LANDAU-LIFSHITZ FORMALISM 

In this section, we briefly review the Landau-Lifshitz 
formulation of gravity and the statement of 4-momentum 
conservation within this theory. Landau and Lifshitz, 
in their Classical Theory of Fields (hereafter referred 
to as LL), reformulated general relativity as a nonlin- 
ear field theory in flat spacetime [49]. (Chap. 20 of 
MTW [50 and a paper by Babak and Grishchuk [ST are 
also helpful sources that describe the formalism.) Lan- 
dau and Lifshitz develop their formalism by first laying 
down arbitrary asymptotically Lorentz coordinates on a 
given curved (but asymptotically- flat) spacetime. They 
use these coordinates to map the curved (i.e. physi- 
cal) spacetime onto an auxiliary flat spacetime (AFS) by 



enforcing that the coordinates on the AFS are globally 
Lorentz. The auxiliary flat metric takes the Minkowski 
form, rj^j, = diag(-l, 1, 1, 1). 

In this formulation, gravity is described by the physical 
metric density 







99'^ 



(2) 



where g is the determinant of the covariant components of 
the physical metric, and g^^ are the contravariant com- 
ponents of the physical metric. When one defines the 
superpotential 



(3) 



the Einstein field equations take the field-theory-in-flat- 
spacetime form 



(4) 



Here r^^ := {-g){T^'" ^t^'^) is the total effective stress- 
energy tensor, indices after the comma denote partial 
derivatives or, equivalent ly, covariant derivatives with 
respect to the flat auxiliary metric), and the Landau- 
Lifshitz pseudotensor tj^^ (a real tensor in the auxiliary 
flat spacetime) is given by Eq. (100.7) of LL [49 or equiv- 
alently Eq. (20.22) of MTW [50j: 



,A0 ^,^Ji 



- ^"^^M-/^p0^^A -/^^M-0"^p0^^A 

+ i (2g'^^g^^ - g'^^g^^) 

X {2gj,pg^r-gpaguT)Q''^,xQ^'',^i (5) 

Due to the symmetries of the superpotential — they are 
the same as those of the Riemann tensor — the field equa- 
tions Q imply the differential conservation law for 4- 
momentum 



T^ru = . 



(6) 



Eq. (|6| is equivalent to T^^ — 0, where the semicolon 
denotes a covariant derivative with respect to the physi- 
cal metric. 

In both LL and MTW, it is shown that the total 
4-momentum of any isolated system (measured in the 
asymptotically flat region far from the system) is 



(7) 



where dT^j is the surface-area element of the flat auxiliary 
metric, and S is an arbitrarily large surface surrounding 
the system. This total 4-momentum satisfies the usual 
conservation law 



dt 



(8) 
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FIG. 4: The regions of space around and inside a binary- 
black- hole system. 



See the end of Section III of [5^ for a brief proof of why 
this holds for black holes. 

Because this paper focuses on BBHs, we will make a 
few further definitions that will be used frequently in our 
study. First, we label the two^ black holes in the binary 
(and the regions of space within their horizons) by A and 
5, and denote their surfaces (sometimes the hole's event 
horizon and other times the apparent horizon) by dA and 
OB, as shown in Fig. [4] We let £ stand for the region out- 
side both bodies but inside the arbitrarily large surface S 
where the system's total momentum is computed (in our 
case, this is taken to be a fixed coordinate sphere inside 
the outer boundary of the numerical-relativity computa- 
tional grid). 

With the aid of Gauss's theorem and the Einstein field 
equations Q, one can reexpress Eq. ^ for the binary's 
total 4-momentum as a sum over contributions from each 
of the bodies and from the gravitational field in the region 
£ outside them: 



Here 



(9a) 



(9b) 



is the 4-momentum of body A (an equivalent expression 
holds for body 5), and 



- f 

Pfleld 



(9c) 



is the gravitational field's 4-momentum in the exterior of 
the black holes. We define an effective velocity of black 
hole A (with similar expressions holding for hole B) by 



Pa 



(10) 



^ After the holes merge, there is only one horizon, which we label 
dC. Equations ([8}-(|10| hold after removing terms with subscript 



In analogy to Eq. ([8| for the rate of change of the bi- 
nary's total 4-momentum, one can write the correspond- 
ing equation for the rate of change of the 4-momentum 
of body A: 



dt 



(r 



_/iO 



dA 



k • 



(11) 



Equation (11) describes the flow of field 4-momentum 



into and out of body A (the second term comes from the 
motion of the boundary of body A with local coordinate 
velocity v\).^ 

We will use Eqs. ([8|-([l0| as the basis for our study 
of momentum flow in black- hole binaries. The actual 
values of the body and field 4-momenta, computed in 
the above ways, will depend on the arbitrary mapping 
between the physical spacetime and the AFS; this is the 
gauge-dependence that will be discussed in Sec. IV B 



III. SIMULATIONS OF HEAD-ON BBH 
COLLISIONS WITH ANTI-ALIGNED SPINS 

In order to investigate the gauge dependence of our re- 
sults, we compare simulations of the same physical sys- 
tem using two separate methods that employ different 
choices of coordinates. One method is a pseudospectral 
excision scheme based on generalized harmonic coordi- 
nates; the other is a finite-difference moving-puncture 
scheme that uses l+log slicing and a gamma-driver shift 
condition (henceforth referred to as "moving puncture 
gauge"; for details see Appendix B2). The coordinates 



used in the two methods differ both for the initial data 
and during the evolution. In this section we summarize 
the construction of initial data and the evolution scheme 
for both methods, and we present convergence tests and 
estimate numerical uncertainties. Further details about 
our numerical methods are are given in Appendices [A] 
and El 



A. Pseudospectral 



Quasiequilibrium excision data 



The evolutions described in Sec. Ill A 2 begin with 
quasiequilibrium excision data constructed using the 
method of Ref. [52 . This method requires the arbitrary 
choice of a conformal three- metric; we choose this metric 



^ In the case that the body's event horizon is stationary (i.e. suf- 
ficiently far from merger), 



dx\ ^y^/dt^ the center of mass 



B and then substituting A 



C. 



velocity of body A. However, if the body's event horizon is dy- 
namical (i.e. during the merger phase), then v\ is the local 
coordinate velocity of the event horizon surface, = dx^^/dt. 
See Sec. |IV A 2| for a discussion of the dynamics of the event 
horizon. 
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SI 


3.902 


0.4986 


0.5162 


±0.5000 


PI 


4.211 


0.4970 


0.5146 


±0.5000 


P2 


8.368 


0.4802 


0.5072 


±0.5091 


HI 


14.864 


0.4870 


0.5042 


±0.5000 



TABLE I: Parameters of the initial data configurations stud- 
ied in this work. Model SI (see Sec. HIAl) gives the pa- 
rameters used to construct a set of Superposed-Kerr-Schild 
quasiequilibrium excision initial data. Model HI (see Ap- 
pendix |d]) gives the parameters for the larger separation 
Superposed-Harmonic-Kerr initial data set. Both SI and HI 
were used in spectral evolutions. PI and P2 provide the 
Bowen-York parameters for the two systems evolved with the 
moving puncture method. The holes are initially separated by 
a coordinate distance d = 2xo and are located at coordinates 
{x,y,z) = (±xo,0,0). For clarity, only 4 significant figures 
are shown. 



to be flat almost everywhere but curved (such that the 
metric is nearly that of a single Kerr-Schild hole) near 
the horizons. 

Our initial data method also requires us to choose an 
outer boundary condition on a shift vector for a gen- 
eral binary that is orbiting and inspiraling, we use^ 



f3' = {noXry^dor'^V^, r 



00, 



(12) 



where is the angular velocity, dor* is the initial radial 
velocity, and Vq is a translational velocity. Note that 
Eq (12) is different from the choice made in Ref. [52] . 



In this paper we confine our focus to collisions that are 
head-on, which we define as l^o = ^0 = 0- However, Vq 
must be nonzero to make the total linear momentum of 
the initial data vanish. 

Table [l| summarizes the initial data used in this pa- 
per. The Arnowitt-Deser-Misner (ADM) mass Madm 
(Eq. (11.2.14) in Ref. [53 ; see also [5l [55]), the irre- 
ducible mass Mirr and Christodoulou mass Mchr of one 
of the holes are listed, where Mchr is related to Mi^r and 
the spin of the hole Sz by 



4M 



2 • 

irr 



(13) 



Table |l] also shows the dimensionless spin Sz/Mq^^^; by 
definition, this measure of the spin lies in the interval 
-1 < 5,/M2,^ < 1. 

For set SI listed in Table El Vq is adjusted so that the 
initial effective velocity of the entire spacetime vl^^ := 
pI^^/p^^^ is smaller than 0.1 km/s, which is approximately 
the size of our numerical truncation error (cf. Fig. H): 



The shift vector used here and in Appendix [a] for the con- 
struction of initial data is not the same as the shift vector used 
during our evolutions. Except for Sec. |III ATj and Appendix [A| 
we always use to refer to the shift during the evolution. 



(Ktl , \vlt\ , Ktl) = (4 X 10-^ 5 X 10-2, 2 X 10-3) km/s 
at time t = 0. 

The construction of initial data is described in more 
detail in Appendix [A] 



2. Generalized harmonic evolutions 

We evolve the quasiequilibrium excision data described 
in Sec. |III A 1| pseudospectrally, using generalized har- 
monic gauge [56, 57, 58, 59 , for which the coordinates 



satisfy the gauge condition 



(14) 



where Hj^ is a function of the coordinates and the space- 
time metric. In this subsection, we summarize the com- 
putational grid used for our spectral evolutions, and we 
briefly discuss our numerical accuracy. Details of our 
pseudospectral evolutions are given in Appendix |B 1[ 

Our computational grid covers only the exterior re- 
gions of the black holes ("black hole excision"): there 
is an artificial inner boundary just inside each appar- 
ent horizon where no boundary condition is needed be- 
cause of causality. The grid extends to a large radius 
'max ^ 400Madm- a set of overlapping subdomains of 
different shapes (spherical shells near each hole and far 
away; cylinders elsewhere) covers the entire space be- 
tween the excision boundaries and r = rmax- 

Because different subdomains have different shapes 
and the grid points are not distributed uniformly, we de- 
scribe the resolution of our grid in terms of the total 
number of grid points summed over all subdomains. We 
label our resolutions NO, Nl, and N2, corresponding to 
approximately 55^, 67^, and 79^ grid points, respectively. 
After merger, we regrid onto a new computational do- 
main that has only a single excised region (just inside the 
newly-formed apparent horizon that encompasses both 
holes). This new grid has a different resolution (and a 
different decomposition into subdomains) from the old 
grid. We label the resolution of the post-merger grid by 
A, and C, corresponding to approximately 63^, 75^, 
and 87^ gridpoints, respectively. We label the entire run 
using the notation ^Nx.y\ where the characters before 
and after the decimal point denote the pre-merger and 
post-merger resolution for that run. Thus, for example, 
'N2.5' denotes a run with approximately 67^ grid points 
before merger, and 75^ grid points afterward. On the out- 
ermost portion of the grid (farther than ~ 200Madm), 
we use a coarser numerical resolution than we do else- 
where. (We only measure the gravitational wave flux, 
linear momentum, etc., at radii of r < 160Madm-) 

To demonstrate the convergence of our evolutions, we 
plot the constraint violation in Fig. [5] for several resolu- 
tions. The quantity plotted is the norm of all the con- 
straints of the generalized harmonic system, normalized 
by the norm of the spatial gradients of all the dynam- 
ical fields, as defined by by Eq. (71) of Ref. [59|. The 
left portion of the plot depicts the constraint violation 
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FIG. 5: Constraint violation at different numerical resolutions 
for the pseudospectral evolutions SI. The common apparent 
horizon forms at time t — 34.73Madm. Labels of the form 
Nx.y indicate the grid resolution, where the pre-merger res- 
olution is labeled (from coarse to fine) by x = 0,1,2 and 
the post-merger resolution is labeled hy y = A, C. The 
constraints decrease exponentially with higher resolution; the 
convergence rate is smaller near merger. 



during the plunge, the right third of the plot shows the 
constraint violation during the ringdown, and the middle 
panel shows the constraints shortly before and shortly af- 
ter the common apparent horizon forms. Throughout the 
evolution, we generally observe exponential convergence, 
although the convergence rate is smaller near merger. 
After merger, there are two sources of constraint vio- 
lations: those generated by numerical truncation error 
after merger (these depend on the resolution of the post- 
merger grid) and those generated by numerical trunca- 
tion error before merger and are still present in the so- 
lution (these depend on the resolution of the pre-merger 
grid). We see from Fig. [s] that the constraint violations 
after merger are dominated by the former source. Also, 
at about t = 200Madm, the constraint violation increases 
noticeably (but is still convergent); at this time, the out- 
going gravitational waves have reached the coarser, out- 
ermost region of the grid. 

Finally, in Fig. [6) we demonstrate the accuracy of the 
recoil velocity 'Ukick = 22 km/s inferred from the gravita- 
tional wave signal ^4, which asymptotically is related to 
the gravitational wave amplitudes /i+ and hx hy 



^4 



(15) 



We extract the spin-weighted spherical harmonic coef- 
ficients of ^4(t) from the simulation as described in 
Ref. [4 , and we integrate these coefficients over time to 
obtain /i^^(t), which are the spin- weighted spherical har- 
monic coefficients of h = h-^ — ihx- For each (^, m), the 
integration constant is chosen so that the average value 
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FIG. 6: Recoil velocity for initial data set SI inferred from the 
gravitational wave signal ^4 extracted at rextr = 160Madm 
at the highest resolution (upper panel). Differences between 
several coarser resolutions and the highest resolution result 
are plotted in the lower panel. 



of h^^{t) is zero. The h^^{t) are then used to compute 
the 4-momentum flux of the gravitational waves from 
Eqs. (3.14)-(3.19) of Ref. [60]. Integrating this flux over 
time yields the total radiated energy- momentum, p^g^^. 
The recoil velocity can then be computed from energy- 
momentum conservation: '^^^^j^ = — pJad/^finab where 
^finai •= ^ADM — ^^rad and £^rad IS the energy radiated 
to infinity. For set SI, we obtain a radiated energy of 
^rad/MADM = (5.6840±0.0008) X 10-"^, where the quoted 
error includes truncation error and uncertainty from ex- 
trapolation to infinite radius (as discussed below). The 
top panel of Fig. |6] shows the recoil velocity as a function 
of time for our highest resolution simulation, while the 
lower panel shows differences between the highest reso- 
lution (N2.C) and lower resolutions. From these differ- 
ences, we estimate a numerical uncertainty for the final 
recoil velocity of 5 x 10~^ km/s for Nl.B and 2 x 10~^ 
km/s for N2.5. 

This numerical uncertainty includes only the effects 
of numerical truncation error; however, there are other 
potential sources of uncertainty in the simulations that 
must also be considered. The first is the spurious "junk" 
gravitational radiation that arises because the initial data 
do not describe a perfect equilibrium situation. This ra- 
diation is not astrophysically realistic, but by carrying a 
small amount of energy-momentum that contributes to 
the measured p^^^ at large distances, the spurious radia- 
tion does affect our determination of the final recoil veloc- 



ity. In our investigation of momentum flow (Sec. IV), we 
do not correct for the initial data's failure to be in equi- 
librium; here we estimate the contribution of the result- 
ing spurious radiation to the final recoil velocity. First, 
we note that for head-on collisions, the physical gravi- 
tational waves are emitted predominantly after merger. 
Therefore, we estimate the influence of the spurious ra- 
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diation by examining the accumulated recoil velocity at 
time t = At + r, where r is the radius of the extraction 
surface and At is a cutoff time. Because the holes merge 
so quickly (because they begin at so small an initial sep- 
aration), the spurious and physical contributions to the 
recoil are not clearly distinguishable in Fig. [6| Vary- 
ing At between 31.1Madm and 38.3Madm (the common 
event and apparent horizons form at t = 31.1Madm and 
t = 34.7Madm, respectively), we estimate that the spu- 
rious radiation contributes approximately 1 km/s (about 
5%) to the recoil velocity — a much larger uncertainty 
than the truncation error. The same variation of At im- 
plies that the spurious radiation contributes about 10% 
of the total radiated energy £^rad-) 

Another potential source of uncertainty in v^-^j^ arises 
from where on the grid we measure the gravitational ra- 
diation. In particular, the quantity ^4 in Eq. ( 15 ) should 



ideally be measured at future null infinity. Instead, we 
measure ^4 on a set of coordinate spheres at fixed radii, 
compute vl^^y^ on each of these spheres, and extrapolate 
the final equilibrium value of v^^^y. to infinite radius. The 

measured from ^4 
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show v^.^^ 



dotted curves on Fig. 
at several radii, and theblack cross shows the final value 
of v^^^y. extrapolated to infinity. We estimate our uncer- 
tainty in the extrapolated value by comparing polynomial 
extrapolation of orders 1,2, and 3; we find an uncertainty 
of 3 X 10~^ km/s for the quadratic fit. Note that if we 
had not extrapolated to infinity, but had instead simply 
used the value of v^^^^. at our largest extraction sphere 
(r = IGOMadm), we would have made an error of 0.85 
km/s, which is much larger than the uncertainty from 
numerical truncation error. Finally, we mention that our 
computation of ^4 is not strictly gauge invariant unless 
^4 is evaluated at future null infinity. As long as gauge 
effects in ^4 fall off faster than 1/r as expected, extrap- 
olation of v^.^j^ to infinity should eliminate this source of 
uncertainty. 



B. Moving puncture 

1. Bow en- York puncture data 

In order to address the importance of gauge depen- 
dence for our calculations using the Landau-Lifshitz for- 
malism, we also simulate BBH mergers using the so- 
called moving puncture method, which employs the co- 
variant form of "1+log" slicing [21 [61] for the lapse func- 
tion a and a "Gamma-driver" condition (based on the 
original "Gamma-freezing" condition introduced in [62j) 
for the shift vector. The precise evolution equations for 
the gauge variables as well as further technical details of 
our puncture simulations are given in Appendix |B 2| 

Our simulations start with puncture initial data [63] 
provided in our case by the spectral solver of Ref. [64|. 
The initial data are fully specified in terms of the initial 
spin Si linear momentum Pi 2 and initial coordinate 
position xi^2 as well as the bare mass parameters mi ^2 of 




FIG. 7: Gravitational recoil for model PI as estimated 
from the gravitational wave signal ^4 extracted at rex = 
73.5 Madm using the highest resolution (upper panel). Dif- 
ferences in the recoil obtained at coarse, medium and fine res- 
olution rescaled for second order convergence (lower panel). 



either hole [65^. The corresponding nonvanishing values 
for the two puncture models considered in this work are 
given m Table |T| There we also list the total black-hole 
mass Mchr and normalize all quantities using the total 
ADM mass Madm- The main difference between the two 
configurations is the initial separation of the holes. The 
lapse and shift are initialized as a = 7"^/^ and = 0, 
where 7 is the determinant of the physical three- metric. 



2. Moving puncture evolutions 

The evolution of the puncture initial data is performed 
using sixth order spatial discretization of the BSSN equa- 
tions combined with a fourth order Runge-Kutta time in- 
tegration. Mesh refinement of Berger-Oliger [66 type is 
implemented using Schnetter's Carpet package f67l[68]. 
The prolongation operator is of fifth order in space and 
quadratic in time. Outgoing radiation boundary con- 
ditions are implemented using second-order accurate ad- 
vection derivatives (see, for example. Sec. VI in Ref. [69]). 

Using the notation of Sec. II E of Ref. [70] the grid 
setup in units of Madm for these evolutions is given by 
(rounded to 3 significant digits) 

{(202,101,58.8,25.2,12.6) x (3.15,1.58,0.788), h}, 
{(201,100,58.5,25.1) x (6.27,3.13,1.57,0.784), h}, 

respectively. Here h denotes the resolution on the in- 
nermost refinement level. For model PI we perform a 
convergence analysis by setting h to he = Madm/49.5, 
hm = Madm /5 7.1 and hf = Madm/64.7, respectively, 
for coarse, medium and fine resolution. Model P2 is 
evolved using h = Madm/49.8. 

Before we discuss the physical results from the punc- 
ture simulations, we estimate the numerical errors due to 
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discretization, finite extraction radius and the presence 
of unphysical gravitational radiation in the initial data. 

In order to study the dependence of the results on res- 
olution, we have evolved model PI of Table [l| using dif- 
ferent resolutions and hf on the finest level and 
correspondingly larger grid spacings by a factor of two on 
each consecutive level. The kick velocity from the high 
resolution simulation, as inferred from the gravitational 
radiation fiux at rex = 73.5 Madm, is shown in the upper 
panel of Fig. [7| The bottom panel shows the differences 
between the velocities obtained at the different resolu- 
tions scaled for second order convergence using a factor 
Q2 = 1.49. By using Richardson extrapolation we es- 
timate the error in the final kick for the fine resolution 
run to be 1 km/s or 5 %. We similarly find overall sec- 
ond order convergence for the velocity derived from the 
components of the Landau-Lifshitz tensor as integrated 
over the apparent horizon. The error in that quantity 
barely varies throughout the entire simulation and stays 
at a level just below A'Ull ^ 50 km/s and 60 km/s for 
fine and coarse resolution respectively. 

The gravitational wave signal is further affected by the 
use of finite extraction radius and linear momentum con- 
tained in the spurious initial radiation. We estimate the 
uncertainty due to the finite extraction radius by fitting 
the final kick velocity obtained for the medium resolution 
simulation of model PI at radii rex = 3 1.5... 94. 5 Madm 
in steps of 10.5 Madm- The resulting final kick veloci- 
ties are well approximated by a polynomial of the form 

+<^i/^ex + <^2/^ex- ^^r = 73.5 M we thus obtain an 
uncertainty of 0.4 km/s corresponding to a relative error 
of 2.2 %. 

Finally we take into account contributions from the 
spurious initial radiation by discarding the wave signal up 
to t— rex = At. For model PI it is not entirely clear where 
exactly the spurious wave signal stops and the physical 
signal starts. By varying At from 30 to 45 Madm we 
obtain an additional error of about ±1 km/s. For model 
P2 no such problem arises because of the smaller am- 
plitude of the spurious radiation and because the longer 
pre-merger time enables the junk radiation to escape the 
system long before the merger happens. We estimate the 
resulting total uncertainty by summing the squares of the 
individual errors and obtain 7.5 % and 5.5 % for models 
PI and P2, respectively. 

Using these uncertainties, the gravitational wave emis- 
sion for model PI results in a total radiated energy of 
^rad/MADM = (0.042 ± 0.008) % and a recoil veloc- 
ity ^kick = (20.3 ± 1.5) km/s. For model P2 the re- 
sult is ^rad/MADM = (0.0555 ± 0.0023) % and ^kick = 
(19.7 ± 1.1) km/s. 



IV. MOMENTUM FLOW 



In this section, we turn to the momentum fiow during 
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FIG. 8: The velocity of the individual and merged black 
holes. The Landau-Lifshitz velocity v^j^ := Pll/^^ll? where 
is ^1^^ Landau-Lifshitz 4-momentum enclosed, is measured 
on the individual and common apparent horizons (labeled AH 
and AHC, respectively) and also on the event horizon (labeled 
EH). For comparison, the coordinate velocities t'coord 
apparent horizons are also shown. The data shown are from 
the high-resolution evolution N2.C. 



merger, and ringdown during a pseudospectral evolution 
of initial data set SI (Table |l|, focusing on the momen- 
tum density and the inferred Landau-Lifshitz velocity 
along and opposite the frame-dragging direction (which 
in this paper a re cho sen to be the direction, respec- 
tively). In Sec. IV we look at the momentum fiow in 
a moving-puncture simulation with similar initial data, 
and by comparing the puncture and spectral simulations, 
we investigate the infiuence of the choice of gauge on our 



results. Then, in Sec. IV C we compare the momentum 
density and velocity of the holes with post-Newtonian 
predictions. 



the evolutions described in Sec. Ill First, in Sec. IV A 



we measure the momentum of the holes during plunge. 



A. Pseudospectral results 



Throughout the pseudospectral evolutions summarized 
in Sec. |III A2[ we measure the 4-momentum density by 
explicitly computing the Landau-Lifshitz pseudotensor 
[Eq. ([5])]. Because our evolution variables are essentially 
the spacetime metric g^jy and its first derivative g^u,p^ 
we are able to compute the momentum density with- 
out taking any additional numerical derivatives. Besides 
measuring the momentum density, we also measure the 
4-momentum [Eq. (9b)] enclosed by i) the apparent 
horizons, ii) the event horizon, and iii) several spheres of 
large radius. From the enclosed momentum, we evaluate 
the effective velocity Vj^^ [Eq. (Il0|]. 



10 




10"-^ 
10' 



INO.B - N2.CI 
INl.B -N2.CI 
IN2.A - N2.CI 
IN2.B - N2.CI 



10 20 



30 34 36 
Time / M 



100 200 300 400 



10^ 
10^ 
^10^ 

10-^ 

10"' 
10"' 
10"' 
10"' 



10" 



ADM 



FIG. 9: Convergence of v^j^ with resolution. Specifically, 
differences between v^-^ at the highest resolution N2.C and at 
various lower resolutions are shown. Labels of the form Nx.y 
indicate the grid resolution, where the pre-merger resolution 
is labeled (from coarse to fine) by x = 0, 1,2 and the post- 
merger resolution is labeled by y = A, B,C. The difference 
between the second-highest and highest resolution is below 
O.lkm/s except near merger, when it grows as large as 1 km/s. 



1. Apparent horizons 

The effective velocities of the apparent horizons are 
shown in Fig.Js] (dashed curves). To demonstrate con- 
vergence, Fig. [9[shows the differences between apparent- 
horizon effective velocities computed at different reso- 
lutions. During the plunge, the difference between the 
medium and fine resolution is less than 0.1 km/s until 
shortly before merger, when it reaches a few tenths of a 
km/s. Shortly after merger, the difference between the 
highest and medium continuation resolutions between 
N2.B and N2.C falls from about 1 km/s to about 0.1 
km/s. 

For comparison. Fig. [S] also shows the apparent hori- 
zons' coordinate velocities (dotted curves); the coordi- 
nate and effective velocities agree qualitatively during 
the plunge and quantitatively during the merger. Also, 
Fig. [8] shows that the effective velocities of individual ap- 
parent horizons and the the event horizon agree well until 
shortly before merger, when the event horizon's velocity 
smoothly transitio ns to ag ree with the common apparent 
horizon's (cf. Sec. IV A 2| below). 

Because of frame-dragging, during the plunge the in- 
dividual apparent horizons accelerate in the downward 
{—y) direction, eventually reaching velocities of thou- 
sands of km/s. But when the common apparent horizon 
appears, its velocity is much closer to zero and quickly 
changes sign, eventually reaching speeds of about 1000 
km/s in the -\-y direction (i.e., in the direction opposite 
the frame-dragging direction). Then, as the common 
horizon rings down, the velocity relaxes to a final kick 



velocity of about 20 km/s in the -\-y direction. 

After merger, why have the horizon velocities suddenly 
changed from thousands of km/s in the frame- dragging 
direction to over a thousand km/s in the opposite direc- 
tion? The answer can be seen in Fig. [lOj which plots con- 
tours of constant y-momentum density at several times. 
At t = 0, the momentum density has an irregular shape, 
because the initial data is initially not in equilibrium. 
By time t = 26.92Madm7 the momentum density has 
relaxed. When the common apparent horizon forms (at 
time t = 34.73Madm), it encloses not only the momen- 
tum of the individual apparent horizons but also the mo- 
mentum in the gravitational field between the holes. 

It turns out that the net momentum outside the indi- 
vidual horizon but inside the common horizon points in 
the -\-y direction; as the common horizon expands, it ab- 
sorbs more and more of this upward momentum. Fig. [iT 
compares the common apparent horizon's effective ve- 
locity to its area and shape; the latter is indicated by 
the pointwise maximum and minimum of the horizon's 
intrinsic scalar curvature. During the first half-period 
of oscillation (to the left of the leftmost dashed vertical 
line), the common horizon expands (as seen by its in- 
creasing area); as it expands, the upward-pointing linear 
momentum it encloses causes v^^ to increase. After the 
first half-period, the horizon shape is maximally oblate 
(cf. panel B on the right side of of Fig. 11), and v^l 
at its maximum value of about 1000 km/s7^ 

After another half-period of oscillation, the apparent 
horizon becomes prolate and encloses enough downward- 
pointing momentum that vf^ has decreased to only about 
+200 km/s. After one additional full period, the effec- 
tive velocity has fallen to nearly zero. As the horizon is 
ringing down, the momentum density in the surrounding 
gravitational field also oscillates: the final four panels 



in Fig. 10 show how the momentum density relaxes to 
a final state as the horizon relaxes to that of a boosted 
Schwarzschild black hole. 

As the horizon rings down, gravitational waves are 
emitted, and these waves carry off a small amount of 
linear momentum. The net radiated momentum is only 
a small fraction of the momenta of the individual holes 
at the time of merger: the final effective velocity of the 
merged hole is about 20 km/s in the upward-pointing di- 
rection, or about 1% of the individual holes' downward 
velocity just before merger. 

Various measures of the final velocity of the merged 
hole are shown in Fig. 
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The kick velocity '^^^^k' 
which is inferred from the outgoing gravitational waves, 
is measured on four coordinate spheres (with radii R of 
IOOMadm, 120Madm, MOMadm, and 160Madm); the 
effective velocity is measured on the same coordinate 
spheres. We find that the effective velocity vf^^ has 
no significant dependence on the radius of the extrac- 
tion surface at late times, while v^^^^ does. The depen- 
dence of v^.^j^ on the extraction radius is expected, since 
our method of extracting ^4 at finite radius has gauge- 
dependent contributions that vanish as ^ 00. When 



11 



Time = 0. Madm 



Time = 26.92 Madm 



Time = 34.73 Madm 






Time = 39.41 Madm 



Time = 4 1 .75 Madm 



Time = 47.6 Madm 






Time = 351.15 Madm 




FIG. 10: Contour plots of the y (up-down) component of the momentum density, which points along or opposite of the holes' 
motion due to frame dragging. Adjacent contours correspond to a factor of 10 difference in the magnitude of the momentum 
density. Contours of positive y momentum density are shown as solid red lines, while contours of negative y momentum density 
are shown as dashed blue lines. The region containing positive y momentum density is shaded grey. The regions inside the 
apparent horizons are shaded black, except for the upper right panel, where the region inside the individual horizons is shaded 
black, while the common apparent horizon is indicated by a thick black line. The data shown are from the high-resolution 
evolution N2.C 



^kick 



is extrapolated to infinite radius^, however, it does 



^ To extrapolate, we fit the velocities v^^^^ at the final time to a 
function of radius R of the form ao -\- ai/R-\- a^ij . 



agree well (within 0.2 km/s) with v\^^. Also, the effec- 
tive velocity v^l calculated on the horizon also agrees 
fairly well (within about 0.5 km/s) with v\^-^ measured 
on distant spheres. 



12 



10^ 



0.6 

u 

S 0.4 

2 



o 

^ 1.9 



1.8 
1.7 



- ; / 




1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 
' ^ 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 i_ 
- 


- 7 








- 


>^ \ B 

- ; \ 


c 


D 


E 




F 

- M^J mm(R)- 

— — — . — 










2 

-- M^^ max(RT 


1 :| 1 1 1 1 1 


1 1 1 1 1 1 


1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


I'l 1 1 1 1 1 


1 1 1 1 1 1 


1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



A 

O 




0.490 0.495 0.500 0.505 0.510 



35 40 45 50 55 60 65 70 75 80 



Dimensionless scalar curvature 
of common apparent horizon (N2.C) 



FIG. 11: Left: A comparison of the common apparent horizon's effective velocity and the horizon's shape and area. The top 
panel shows the horizon's effective velocity v^-^. The middle panel shows the pointwise minimum and maximum of the horizon's 
dimensionless intrinsic scalar curvature; both MQi^j,mm{R) and Mq^^^, msix{R) relax to the Schwarzschild value of 1/2 as the 
horizon rings down. (The first four local minima of M^^r niin(i?) are indicated by vertical dashed lines.) The bottom panel 
shows the area A of the common apparent horizon normalized by the total area of the individual horizons at t = 0. The data 
shown are from the high-resolution evolution N2.C. Right: The dimensionless intrinsic scalar curvature Mqy^^R of the common 
apparent horizon at the times labeled A-F in the left panel. The horizon begins peanut-shaped, then rings down, eventually 
settling down to a sphere with a constant curvature M^^^^R = 0.5. 
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FIG. 12: A comparison of various measures of the merged 
hole's final velocity. The kick velocity inferred from the 
gravitational-wave flux (thin dashed lines) and the Landau- 
Lifshitz effective velocities v\^^ (thin solid lines) are measured 
on spheres of radius IOOMadm, 120Madm, 140Madm, and 
160Madm The value of the kick velocity at the final time is 
extrapolated to r = oo (black cross). The effective velocity 
measured on the common apparent horizon (thick solid line) 
and the coordinate velocity (thick dashed line) are also shown. 
The data shown are from the high-resolution evolution N2.C 
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FIG. 13: The effective velocity v\^^ calculated on the event 
horizon surface, with the specified snapshots in Fig. [14] of 
the event horizon surface marked: a,b, t = 27.7Madm; c, 
t 30.8Madm; d, t 31.6Madm; e, t 35.5Madm; f, 
t = 40.8Madm. 



2. Event horizon 

We would like to compare our quantitative results of 
the effective velocity calculated using the event hori- 
zon surface (Fig. 13) with qualitative observations of 

We find that 



the event horizon's dynamics (Fig. 14). 
the greatest variation in both the event horizon geome- 
try and the value of v\^^ occurs over a period of about 
At 13Madm from t = 28Madm to t = 41Madm. 
At time t = 27.7Madm, the cusps of the event hori- 
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FIG. 14: Snapshots of the event horizons at the times indicated in Fig. [13] a,b, t = 27.7Madm; c, t = 30.8Madm; d, 
t = 31.6Madm; e, t = 35.5Madm; f, ^ = 40.8Madm- All snapshots are looking down the z-axis to the x-y plane, except for 
shot a, which is slightly skewed (slightly rotated about the y axis) to better see the geodesic structure. In shot a, the future 
generators of the horizon are visible as small blue dots. Note how the future generators map out a surface that meets the 
event horizon at the event horizon's cusps; this is where the future generators join the horizon. The data shown are from the 
high-resolution evolution N2.C. 



zon just begin to become noticeable (Figs. 14 a & b). 
One can see in Fig. 13 that this is the time at which 



changes from decreasing to increasing. Shortly af- 
ter^, at t = 31.1Madm, the two separate event horizons 



coalesce into a common event horizon, and the common 
event horizon rapidly expands to form a convex shape 
by t = 35.5Madm (Figs. [l4]d & e). At this time, we 
note that vf^j^ is rapidly increasing (Fig. [13] arrow e); 
this rapid increase corresponds to the quickly expanding 
event horizon surface. 



^ Note that at t = SI .IMadm^ we (smoothly) modify our gauge 
condition [Eq. ( |B11| and the surrounding discussion]. The sepa- 
rate event horizons coalesce at time t = 31.1Madm as well; this 
is a coincidence. 



We interpret this process as the merging black holes 
"swallowing" the gravitational field momentum between 
the holes. The resulting change in I'll divided 
into two distinct portions: i) one that results from the 
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changing event horizon surface in space, i.e. the field mo- 
mentum swahowed by the black holes [mathematically, 
the second term, in Eq. (11)] and ii) a second that re- 
sults from the change of field momentum at the black 
holes' surface, i.e. the field momentum flowing into the 
black holes [mathematically, the first term, in Eq. (11)]. 
While this distinction is clearly coordinate dependent, it 
could, after further investigation, nevertheless provide an 
intriguing and intuitive picture of the near-zone dynam- 
ics of merging black hole binaries. 



B. Moving- puncture results and gauge 

As summarized in Sec. [llj the Landau-Lifshitz formal- 
ism that we have applied to our numerical simulations 
is based on a mapping between the curved spacetime of 
the simulation and an auxiliary fiat spacetime. In the 
asymptotically-fiat region far from the holes, there is a 
preferred way to construct this mapping. Consequently, 
when the surface of integration is a sphere approaching 
infinite radius, Eq. (9b) gives a gauge-invariant measure 



of the system's total 4-momentum (see, e.g.. Sec. 20.3 of 
Ref. [50 ). However, when the surface of integration is in 
the strong-field region of the spacetime (e.g., when the 
surface is a horizon), the 4-momentum enclosed is gauge 
dependent. The momentum density, being given by a 
pseudotensor, is always gauge dependent. 

The gauge-dependence of the effective velocity can be 
investigated at late times — when the spacetime has re- 
laxed to its final, stationary configuration — by comparing 
the velocity obtained on the horizon with gauge- invariant 
measures of the kick velocity (Fig. 12). At the final time 
in our pseudospectral simulation, the effective velocities 
of the apparent and event horizons agree within tenths 
of a km/s with the (extrapolated) kick velocity inferred 
from the gravitational- wave fiux; at late times, the hori- 
zon effective velocities also agree with the effective veloc- 
ity measured on coordinate spheres of large radius. At 
least at late times, then, the effective velocity I'll 
significantly affected by our choice of gauge. 

But how strong is the infiuence of gauge on our results 
in the highly-dynamical portion of the evolution, when 
we have no gauge- invariant measure of momentum or ve- 
locity? To investigate this, we have evolved initial data 
that are physically similar using two manifestly differ- 
ent gauge conditions: i) the generalized-harmonic condi- 
tion used in our spectral evolutions, and ii) the "1+log" 
slicing and "Gamma-driver" shift conditions used in our 
moving-puncture evolutions. 



Figs. [15] and [16] display the velocity obtained from the 
horizon integral of the components of the Landau-Lifshitz 
tensor in the moving-puncture evolutions described in 
Sec. Ill B 2 The most remarkable feature in these plots 





is a large temporary acceleration of the black holes in the 
frame-dragging direction. The magnitude of the velocity 
reaches about 4500 km/s, which is of the order of the 
superkicks first reported in Refs. |30l[32]. In contrast to 



FIG. 15: Velocity perpendicular to the line of sight associ- 
ated with the horizon integrals of the the Landau-Lifshitz 
tensor obtained for model PI. The shaded area represents 
the numerical uncertainty. During the pre-merger phase, the 
velocities of both holes are identical. 
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FIG. 16: Same as Fig.[T5|for model P2 of Table|lj 



those inspiraling configurations, however, the black hole 
motion reverses during the merger and settles down to a 
small value of —30 ± 50 km/s. 

In order to examine to what extent this behavior is 
dependent on specific properties of the puncture evolu- 
tion (such as the particular form of the spurious radia- 
tion, which differs in our spectral and puncture evolu- 
tions), we have performed the following additional simu- 
lations. First, we have changed the gauge parameter 77 in 
Eq. ( |B20 ) to 0.75 and 1.25. We do not observe a signif- 
icant change in the behavior of the effective velocity for 
this modification. 

Second, in order to gain further insight into the depen- 
dence of the effective velocity on the initial separation of 
the black holes, we have increased the initial separation 
of the holes to allow for a longer pre-merger interaction 
phase; We study the evolution of the second model P2 
in Table [l] This simulation has been performed with the 
Lean code as summarized in Sec. [Ill BT] using a resolu- 
tion he = Madm/49.8. The resulting velocity is shown 
in Fig. [16] and represents numerical uncertainties as gray 
shading. The remarkable similarity between the figure 
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and its counterpart Fig. [15] for model PI demonstrates 
that the numerical results are essentially independent of 
the initial separation. 

Comparing Figs. [S] and [T5| the qualitative behavior of 
the apparent horizons' effective velocities agrees. In both 
the spectral and puncture simulations: 

1. during the plunge, the individual apparent horizons 
accelerate to speeds larger than 1000 km/s in the 
frame dragging direction, 

2. when the common horizon forms, its velocity is 
much smaller in magnitude, because the common 
horizon has enclosed momentum pointing opposite 
the frame-dragging direction, and 

3. the velocity relaxes to a value of only tens of km/s 
that (within numerical uncertainty) agrees with the 
kick velocity measured using the gravitational- wave 
flux. 

These results are particularly encouraging because two 
popular gauge choices used in the NR community give re- 
markable overall agreement. While this qualitative agree- 
ment certainly does not constitute a proof of a gauge 
independence of our findings, we feel encouraged in our 
hope that different types of observers might agree on their 
overall perception of the local black-hole dynamics dur- 
ing the collision. Most importantly from a practical point 
of view, it appears possible that such local descriptions 
can be derived from the current generation of BBH codes 
without the different numerical relativity groups having 
to agree upon one and the same gauge choice for com- 
paring their momentum densities and effective velocities. 
Future investigations using a wider class of coordinate 
conditions should further clarify the significance of gauge 
choices in this context. 



C. Comparison with post-Newtonian predictions 

In this section we compare our results to post- 
Newtonian predictions. For each comparison, first the SI 
data set (Table|l| is presented along with post-Newtonian 
predictions of a corresponding initial configuration, then 
the HI data set (Table |l| is presented along with its post- 
Newtonian predictions. The post-Newtonian trajectories 
for spinning point particles were generated by evolving 
the post-Newtonian equations of motion [TlJ [72]. The 
difference between the two data sets are: i) set HI be- 
gins with a larger initial separation than set SI, and ii) 
set HI is evolved in a nearly harmonic gauge. Comparing 
evolutions of data sets SI and HI illustrates how these 
two effects improve the comparisons one can make with 
post-Newtonian predictions. 

The left panels of Figs. p!7}|l9] show the comparison be- 
tween the highest-resolution evolution (N2.C) of initial 
data set SI and several orders of post-Newtonian predic- 
tions. The right panels of Figs. 17-19 show analogous 



Figure [17] shows that the bulk, longitudinal motions 
(i.e., motion in the x direction) agree both qualita- 
tively and quantitatively with post-Newtonian predic- 
tions through most of the plunge (i.e., a few Madm before 
the formation of the common apparent horizon) for both 
data sets. In the left panel of Fig. [iT] we have added 
another 2.5 PN curve that is offset vertically such that 
the 2.5 PN coordinate velocity agrees exactly with the 
numerical effective velocity at t 18.34Madm; this is 
done in order to account for the period of initial relax- 
ation in the SI data set. Quantitative agreement is then 
found between 2.5 PN predictions and both the effec- 
tive and coordinate velocities from t ^ 5Madm through 
t ^ 20Madm- The right panel of Fig. [iT] which has less 
of an initial relaxation due to the increased separation, 
shows excellent agreement between both the effective and 
coordinate velocities and the 2.0 PN and 2.5 PN predic- 
tions. 

For the minor (yet more interesting) transverse mo- 
tion (i.e., the motion along the y direction), we find only 
qualitative agreement between the numerical data and 
post-Newtonian predictions — spin-orbit coupling [more 
specifically, frame-dragging plus spin-curvature coupling, 
see Eq. (5.11) of paper I and discussions thereafter] 
cause the holes to move in the —y direction during the 
plunge, reaching speeds of order 1000 km/s before the 
holes merge. The post-Newtonian expansion scheme we 
adopt (paper I and Refs. [7Tl|72]) uses a harmonic gauge, 
and a physical spin supplementary condition (SSC) of 
S'^^up = 0, where S'^^ is the spin angular momentum 
tensor of the black hole and its four velocity (see e.g.. 
Sec. II B of paper I). 

In this scheme, for the equal-mass-opposite-spin con- 
figuration, up to the leading 1.5 PN order, the coordinate 
y velocity of the point particle representing each hole is 



equal to 3/2 the hole's effective velocity. 



evaluated 



through a surface integral of the post-Newtonian expres- 
sion for the super potential [cf. Eq. (9b)]. Therefore, in 

If 



Figs. 

factor 
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comparisons with an evolution of initial data set HI. 



we rescale the effective velocity ^^^^ 
3/2, which arises from our particular choice of 
SSC and from field momentum distribution in the vicin- 
ity of the holes (see Sees. II B and II C, and Table I of 
paper I for details). 

In Figs. [18] and [19] we compare the post-Newtonian 
point-particle y velocity with the numerical coordinate 
y velocity and 3/2 of the numerical effective y velocity 
^LL- -^^^ comparison to the SI data set, we find qual- 
itative agreement with both the effective and coordinate 
velocities and the post-Newtonian predictions. We think 
this agreement is not better because of the large initial 
relaxations present in the SI data set related to small 
initial separation. However, in the HI comparison, we 
find excellent agreement between the coordinate velocity 
and the 2.5 PN prediction but only qualitative agree- 
ment between the effective velocity and post-Newtonian 
predictions. In these figures, offsets of —433 km/s (for 
SI data) and —38 km/s (for HI data) have been used to 
make 2.5 PN coordinate velocity agree better with nu- 
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FIG. 17: A comparison of numerical and post-Newtonian longitudinal velocities (i.e., j c) versus time. The predicted 
coordinate velocities at several post-Newtonian orders are shown as broken curves. Left: A comparison of SI numerical data 
and post-Newtonian predictions. The numerical and post-Newtonian curves agree qualitatively. When the 2.5 PN curve is 
offset by a certain amount, it agrees quantitatively with the coordinate velocity t'coord ^^^^ the effective velocity v^i^. Right: A 
comparison of HI numerical data and PN predictions. The effective velocity v^j^ (thick black line) closely tracks the coordinate 
velocity t'coord5 both numerical curves also agree well with the 2.0 PN and 2.5 PN curves. 
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FIG. 18: A comparison of numerical and post-Newtonian transverse velocities (i.e., v'^ in km/s) versus time. The left panel 
shows numerical results from simulation SI, while the right panel shows numerical results from simulation HI. The predicted 
coordinate velocity at several post-Newtonian orders are shown as broken curves. The effective velocity is shown in black; it 
has been rescaled by a factor of 3/2 in order to aid comparison with the post-Newtonian point-particle velocities, as discussed 
in the text. 



merical results. Such offsets can be motivated as fol- 
lows. Our numerical initial data were chosen such that 
the initial total momentum of the entire spacetime van- 
ishes. This, in our post-Newtonian scheme, corresponds 
to nonvanishing initial y velocities of (see Table I of paper 

I) 

" 4(ro/MADM)2 ' ^^^^ 

where x is the spin parameter of each hole, and tq their 
initial separation. This corresponds to —616 km/s for the 
SI data, and —42 km/s for HI data. Again, the agree- 
ment is qualitative for SI data, and quantitative for HI 
data. 

One final comparison we make between the HI data 
set and post-Newtonian predictions is the near-field mo- 
mentum density, shown in Fig. |20j The numerical data 



comes from the harmonic evolution HI, while the 1.5 
PN momentum density is computed from Eqs. (A2a)- 
(A2c) in paper I using the numerical hole trajectories. 
The left panels, comparing the initial data to the pre- 
dicted post-Newtonian momentum density, show differ- 
ences which are presumably due to differences in the 
post-Newtonian and numerical initial data, such as the 
numerical initial data being out of equilibrium. The cen- 
ter panels show the momentum densities agree very well 
once enough time has elapsed for the spacetime to re- 
lax and for the spurious radiation to be emitted but be- 
fore the holes have fallen too close together. The right 
panels make a final comparison just before the holes get 
close enough to merge and shows differences appearing 
between the numerical data and the post-Newtonian pre- 
dictions very near the holes — which could be an indica- 
tion of the breakdown of the post-Newtonian approxima- 
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FIG. 19: A comparison of numerical and post-Newtonian velocities. In the figure, v'^ in km/s is plotted against /c. The 
effective velocity v^^ of the highest-resolution (N2.C) evolution of initial data SI (Table [l]) on the left and of the evolution of 
initial data HI (Table[l| on the right are shown as a thick black line. The predicted coordinate velocity at several post-Newtonian 
orders are shown as broken curves. The transverse effective velocities only agree qualitatively with post-Newtonian predictions; 
however, the coordinate velocity agrees very well with post-Newtonian predictions. In the left panel, the coordinate velocity 
has been artificially truncated shortly before merger because at that point we do not have a good measure of the coordinate 
velocity. The effective velocity has been rescaled by a factor of 3/2 to aid comparison with the post-Newtonian point-particle 
velocities, as discussed in the text. 
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FIG. 20: Comparison of numerical (top row) and post-Newtonian (bottom row) y momentum density. The numerical data 
comes from the harmonic evolution HI described in Appendix |d] The 1.5 PN momentum density is computed from Eqs. (A2a) 



(A2c) in paper I using the numerical hole trajectories. As in Fig. 10 contours represent powers of 10 in y momentum density. 
The positive y momentum density contours are shown in red, negative in blue. The region of positive y momentum density is 
shaded grey. In the numerical plots the apparent horizons are shown in black. 
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tion. 

These comparisons with post-Newtonian predictions 
have yielded several interesting results. The primary re- 
sult of these comparisons is the surprisingly good agree- 
ment found between post-Newtonian predictions and the 
coordinate velocities, especially from the harmonic gauge 
evolution. Also, the longitudinal effective and coordi- 
nate velocities track each other; consequently, the lon- 
gitudinal effective velocity agrees with post-Newtonian 
predictions. The transverse effective velocities agree 
qualitatively with the post-Newtonian predictions in the 
sense that they both indicate that the holes accelerate 
in the expected frame-dragging direction to speeds of or- 
der 1000 km/s. Finally, we have also found the qualita- 
tive agreement between harmonic gauge numerical data 
and post-Newtonian extends to the near-zone momentum 
density after the initial data relaxes but before the holes 
have fallen too close together. 



V. CONCLUSION 

With the goal of building up greater physical intuition, 
we have used the Landau-Lifshitz momentum-flow for- 
malism to explore the nonlinear dynamics of fully rel- 
ativistic simulations of a head-on BBH plunge, merger, 
and ringdown. We have defined and computed an ef- 
fective velocity of the black holes in terms of the mo- 
mentum and mass-energy enclosed by their horizons, and 
we have interpreted the holes' transverse motion — which 
reaches speeds of order 1000 km/s — as a result of momen- 
tum flow between the holes and the gravitational field 
of the surrounding spacetime. We have found that the 
merged hole's final effective velocity — about 20 km/s — 
agrees with the recoil velocity implied by the momentum 
carried off by the emitted gravitational waves. 

Our measures of linear momentum and effective ve- 
locity are gauge dependent. Nonetheless, after compar- 
ing simulations of comparable initial data in generalized- 
harmonic and moving-puncture gauges, we have observed 
remarkably weak gauge dependence for the generalized- 
harmonic and moving-puncture evolutions discussed in 
this paper. Additionally, we have found surprisingly good 
agreement between the holes' effective and coordinate ve- 
locities, and at late times, the holes' final effective veloc- 
ities and gauge- invariant measures of the kick velocity 
agree. 

These results motivate future explorations of momen- 
tum flow in fully-relativistic numerical simulations that 
are more astrophysically realistic. We are particularly 
eager to investigate simulations of superkick BBH merg- 
ers (the inspiral of a superkick configuration was consid- 
ered using the post-Newtonian approximation in paper 
I). Other future work includes studies of the linear and 
angular momentum flow in inspiraling (rather than head- 
on) mergers as well as mergers with larger spins. 
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APPENDIX A: SUPERPOSED-KERR-SCHILD 
(SKS) INITIAL DATA 

The initial data for the pseudospectral simulations pre- 
sented in this paper was constructed using the methods 
described in Ref. [52 . In this appendix, we describe in 
more detail these initial data (which we summarize in 
SeclIIIAl). 



The usual 3+1 decomposition splits the spacetime met- 
ric g^i, into a spatial metric 7^^, lapse a, and shift i.e. 



jij{dx' ^ f3'dt){dx^ + 



P^dt). 



(Al) 



On the initial spatial slice (at time t = 0), the initial 
data must specify the spatial metric jij and the extrinsic 
curvature Kij^ which is related to the time derivative of 



the spatial metric by 



-2aKij 



(A2) 



We use the quasiequilibrium formalism [73j [TH [75j [76l 
.77j . in which 7^^ and Kij are expanded as 



K — A - ■ 



(A3) 



The conformal metric 7^^, the trace of the extrinsic cur- 
vature and their time derivatives can be chosen freely. 
We adopt the quasiequilibrium choices 



dtK 



0, 
= 0. 



(A4) 



The remaining free data are based on a weighted super- 
position of two boosted, spinning Kerr-Schild black holes 
(Eqs. (45)-(46) of Ref. [52,): 



0=1 



a=l 



(A5) 



(A6) 
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Here fij is the metric of fiat space, is the Euchdean 
distance from the center of the apparent horizon of hole 
a, and jfj and Ka are the spatial metric and mean cur- 
vature of a boosted (with velocity t;*), spinning (with 
spin S/M'^) Kerr-Schild black hole centered at the ini- 
tial position of hole a. In this paper we choose i)* = 
(since we seek data describing holes falling head-on from 
rest), M = 0.39Madm, and S/M'^ = 0.5. The Gaussian 
weighting parameter is chosen tohe Wa = d/3^ where d is 
the initial coordinate separation between the two holes; 
note that this choice causes the conformal metric to be 
fiat everywhere except near each hole. The holes are lo- 
cated at coordinates {x^y^z) = {xq = ±(i/2,0,0). 

These free data are then inserted into the ex- 
tended conformal thin sandwich (XCTS) equations (e.g., 
Eqs. (13)-(15) of Ref. [74 )^, which are then solved for 
the conformal factor the lapse a, and the shift 



1 

12 
2 
3 



1 



0, 



r R .1^ 

{aip) 

-i^'idtK-p'dkK), 



0, 



1 + ^x2^4+7^-8^,,-^ 



(A7) 



where the Ricci scalar curvature of the conformal metric 
jij is R, the longitudinal derivative L is defined as 



{lV)i 



^:lij^kV\ (A8) 



and the trace-free part of the extrinsic curvature A^^ sat- 
isfies 



The XCTS equations are solved using a spectral ellip- 
tic solver [78 on a computational domain with i) a very 
large outer boundary (which is chosen to be a coordinate 
sphere with radius lO^M), and ii) with the region inside 
the holes' apparent horizons excised. The excision sur- 
faces S are surfaces of constant Kerr radius rxerr, where 



~. 2 



1. 



^Kerr + /Ma ^Kerr 



(A9) 



The excision surfaces are the apparent horizons of the 
holes; this is enforced by the following boundary condi- 



tion: (Eq. (48) of Ref. [74]): 



(L/3),, 
1 



6 



i^V^^ on S. (AlO) 



Here := is a unit vector normal to the excision 

surface 5, and hij := — SiSj is the induced metric on 
the excision surface. 

On the apparent horizon, the lapse satisfies the bound- 
ary condition 



a^l) = 1 + 



[aa - 1) on 5, 



(All) 



where is the lapse of the Kerr-Schild metric corre- 
sponding to hole a. The shift satisfies 



(5'^ = as^ — oil ^- 



(A12) 



The first term in Eq. (A12) implies that the holes are 



initially at rest, and the second term determines the spin 
of the hole; to make the spin point in the ±z direc- 
tion with magnitude S/Mq^^^ = 0.5 (measured using the 
method described in Appendix A of Ref. [52 ), we choose 
MADM^r = tO- 244146 and = 9^, where is the ro- 
tation vector on the apparent horizon corresponding to 
rotation about the +z axis. 

On the outer boundary S, the spacetime metric is fiat: 



ij) — 1 on S, 
ail) — 1 on S. 



(A13) 
(A14) 



Because the holes are initially at rest in the coordinates, 
they can be given orbital, radial, and translational mo- 
tion by rotation, expansion, and translation of the shift 
on the outer boundary, i.e. 



(3' = {no X rY 



FJ, on B. 



(A15) 



We choose do = and 1] = 0. To make the total momen- 
tum of the initial data vanish, we choose = —0.001444. 

Our initial data are constructed [Eq. ( |A12 )] in a frame 
comoving with the black holes. Thus, an asymptotic ro- 
tation, expansion, and translation in the comoving shift 
cause the holes to initially have radial, angular, or 
translational velocity in the inertial frame. Note that 
the initial data are evolved in inertial, not comoving, co- 
ordinates, so that the shift during the evolution is differ- 
ent from the comoving shift obtained from the XCTS 
equations: the former asymptotically approaches zero, 
not a constant vector Vq. 



APPENDIX B: NUMERICAL METHODS FOR 
EVOLUTIONS 



^ The XCTS equations are also given by Eqs. (37a)-(37d) of 
Ref. 52 , aside from the fohowing typographical error: the sec- 
ond term in square brackets on the right-hand-side of Eq. (37c) 
should read {5 / 12) K'^ i;"^ (not (5/12)X4V^4). 



1. Pseudospectral evolutions 

We evolve the initial data summarized in Sec. IIII A II 
using the Caltech-Cornell pseudospectral code SpEC. 
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This code and the methods it employs are described in 
detail in Refs. [H [79l [80]. Some of these methods have 
been simplified for the head-on problem discussed here, 
and others have been modified to account for a nonzero 
center-of-mass velocity, so we will describe them here. 

We evolve a first-order representation [59] of the gen- 
eralized harmonic system [56[ [57l |58]. We handle the 
singularities by excising the black hole interiors from 
the computational domain. Our outer boundary condi- 
tions [59l [m [82] are designed to prevent the infiux of un- 
physical constraint violations jS^} [84l [85l [86l ISTj [88} [89] 
and undesired incoming gravitational radiation [90l [91] 
while allowing outgoing gravitational radiation to pass 
freely through the boundary. 

We employ the dual-frame method described in 
Ref. [79 : we solve the equations in an "inertial frame" 
that is asymptotically Minkowski, but our domain de- 
composition is fixed in a "comoving frame" that is al- 
lowed to shrink, translate and distort relative to the in- 
ertial frame. The positions of the centers of the black 
holes are fixed in the comoving frame; we account for 
the motion of the holes by dynamically adjusting the co- 
ordinate mapping between the two frames. Note that 
the comoving frame is referenced only internally in the 
code as a means of treating moving holes with a fixed 
domain. Therefore all coordinate quantities (e.g. black 
hole trajectories) mentioned in this paper are inert ial- 
frame values unless explicitly stated otherwise. 

The mapping from comoving to inertial coordinates is 
changed several times during the run. During the plunge 
phase, we denote the mapping by A^p(x*,x'*), where 
primed coordinates denote the comoving frame and un- 
primed coordinates denote the inertial frame. Explicitly, 
Mp{x\ x^'^) is the mapping 



X = F {r\ t) sin 0' cos <p\ 

y = t) sin sin (^' + e-^"/^Ty(t), 

z = F{r' ,t) cos e' cos (j)' , 



where 



F{r\t) :-- 



a(t) + (1 - a{t)) 



p/2 



(Bl) 
(B2) 
(B3) 



(B4) 



Here a{t) and Y{t) are functions of time, {r'^O'^cj)') are 
spherical polar coordinates in the comoving frame cen- 
tered at the origin, and R'q and are constants. For the 
choice Rq = oo and = oo, the mapping is simply an 
overall contraction by a{t) < 1 plus a translation Y{t) in 
the y direction. Choosing Rq equal to the outer bound- 
ary radius i^^ax choosing ~ ^max/^ causes the 
map to approach the identity near the outer boundary; 
this prevents the outer boundary from falling close to 
the strong-field region during merger, and makes it eas- 
ier to keep the outer boundary motion smooth through 
the merger /ringdown transition. The functions a{t) and 
Y{t) are determined by dynamical control systems as 
described in Ref. ^79j. These control systems adjust 



a{t) and Y{t) so that the centers of the apparent hori- 
zons remain stationary in the comoving frame. For the 
evolutions presented here, we use Rq = 532.2Madm = 
l.lR'^g^^ and = 31.21Madm = 4(io, where do is the 
initial separation of the holes. 

The gauge freedom in the generalized harmonic system 
is fixed via a freely specifiable gauge source function Ha 
that satisfies the constraint 







(B5) 



where r^5c are the spacetime Christoffel symbols. To 
choose this gauge source function, we define a new quan- 
tity Ha that transforms like a tensor and agrees with Ha 
in inertial coordinates (i.e. Ha = Ha). Then we choose 
Ha so that the constraint (B5) is satisfied initially, and 
we demand that Ha' is constant in the moving frame. 

Shortly before merger (at time ti = 31.1Madm), we 
make two modifications to our algorithm to reduce nu- 
merical errors and gauge dynamics during merger. First, 
we begin controlling the size of the individual apparent 
horizons so that they remain constant in the comoving 
frame, and therefore they remain close to their respective 
excision boundaries. This is accomplished by changing 
the map between comoving and inertial coordinates as 
follows. We define the map MaUi{x\x''^) for black hole 
1 as 

X = Xahi + ^sin6>' cos^^ (B6) 
y = ^AHi +^sin6>'sin^', (B7) 
z = ^AHi +rcos(9', (B8) 
f := /-e-("'-^o)V-?Ai(t), (B9) 

where {r' ^0' ^(j)') are spherical polar coordinates centered 
at the (fixed) comoving-coordinate location of black hole 
1, which we denote as (^aHi ' ^aHi ' ^aHi)- constant 
^AHi desired average radius (in comoving coor- 

dinates) of black hole 1. Similarly, we define the map 
A1ah2(^%^'*) foi" black hole 2. Then the full map from 
the comoving coordinates x''^ to the inertial coordinates 
X* is given by 

M,n{x\x'') := M^{x\x')MAnAx\x')MAnA^\x'')' 

(BIO) 

The constants ai, cr2, and Tq are chosen to be 
0.780Madm, 0.780Madm, and I.OIMadm, respectively. 
The functions Xiit) and A2(t) are determined by dynam- 
ical control systems that drive the comoving-coordinate 
radius of the apparent horizons towards their desired val- 
ues i^AHi — ^AH2 ~ 1-56Madm Note that in comoving 
coordinates, the shape of the horizons is not necessar- 
ily spherical; only the average radius of the horizons is 
controlled. 

The second change we make at time ti = 31.1Madm 
is to smoothly roll gauge source function Ha to zero by 
adjusting Ha'it) according to 



Ha'{t)=Ha'{tl)e 



-(t-tl)VT^ 



(Bll) 
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where r = 0.5853Madm- This choice makes it easier for 
us to continue the evolution after the common horizon 
has formed, and it also reduces gauge dynamics that oth- 
erwise cause oscillations in the observed Landau-Lifshitz 
velocity 'U^l during the ringdown. 

When the two black holes are sufficiently close to one 
another, a new apparent horizon suddenly appears, en- 
compassing both black holes. At time tm = 34.73Madm 
(which is shortly after the common horizon forms), we in- 
terpolate all variables onto a new computational domain 
that contains only a single excised region, and we choose 
a new comoving coordinate system so that the merged 
(distorted, pulsating) apparent horizon remains spheri- 
cal in the new comoving frame. This is accomplished in 
the same way as described in Section II. D. of [4 , except 
that here the map from the new comoving coordinates to 
the inertial coordinates contains an additional transla- 
tion in the y direction that handles the nonzero velocity 
of the merged black hole. In [4 a third change, namely 
a change of gauge, was necessary to continue the sim- 
ulation af ter m erger. But in the simulations discussed 
here, Eq. ( Bll[ ) has caused Ha to fall to zero by the time 
of merger, and we find it suffices to simply allow Ha to 
remain zero after merger. 

For completeness, we now explicitly describe the map 
from the new comoving coordinates x"'^ to the inertial 
coordinates x*. This map is given by 



X = 

y = 
z — 



r sin 6" cos 

rsmO" sin0'' + e-^"'/^TV(t), 
r cos 6'\ 

l + sin2(7rf/2i?:;,J 



(B12) 

(B13) 
{B14) 



R" 



{1-A{t)) 



ID// D/2 



1 (515) 



E ^Ut)YU0",4>"), (B16) 



i=0 m=-i 



{r"^0"^(j)") are spherical polar coordinates in the new 
comoving coordinate system, i^^ax value of r" at 

the outer boundary, and is a constant chosen to be 
31.21Madm- The function q{r") is given by 



q{r" 



-(»•"- 



(B17) 



where R'xb. is the desired radius of the common apparent 
horizon in comoving coordinates. The function A(t) is 



A(t)=A^^{Ai^A2(t-tm))e 



(Bl^ 



where the constants Aq, Ai, and A2 are chosen so that 
A{t) matches smoothly onto a{t) from Eq. ( |B4| ): A(tm) = 

a(tm)^ A{tm) = Ciitm)^ and A(tm) = o^itm)' The constant 
ta is chosen to be on the order of 5M. The functions Y{t) 
and A^^(t) are determined by dynamical control systems 
that keep the apparent horizon spherical and centered at 
the origin in comoving coordinates; see [4| for details. 



2. Moving- puncture evolutions 

In addition to the spectral evolutions, we have per- 
formed a second set of simulations using the so-called 
moving puncture technique [2] [3] using the Lean code 
[701 [92] . This code is based on the Cactus computa- 
tional toolkit [93 and uses mesh refinement provided by 
the Carpet package [67jj68^. Initial data are provided 
in the form of the TwoPUNCTURES thorn by Ansorg's 
spectral solver [64 and apparent horizons are calculated 
with Thornburg's AHFinderDirect [94. .95^ . 

The most important ingredient in this method for the 
present discussion is the choice of coordinate conditions. 
A detailed study of alternative gauge conditions in the 
context of moving puncture type black-hole evolutions is 
given in Ref. [96]. In particular, they demonstrate how 
the common choice of a second order in time evolution 
equation for the shift vector can be integrated in time 
analytically and thus reduced to a first order equation. 
Various test simulations performed with the Lean code 
confirm their Eq. (26) as the most efficient method to 
evolve the shift vector. In contrast to the shift, moving 
puncture codes show little variation in the evolution of 
the lapse function. Here we follow the most common 
choice so that our gauge conditions are given by 



dtp' 



(B19) 



(B20) 



is the contracted Christoffel symbol of the conformal 
3-metric, K the trace of the extrinsic curvature [see for 
example Eq. (1) of [70 J and 77 a free parameter set to 
1 unless specified otherwise. For further details about 
the moving puncture method and the specific implemen- 
tation in the Lean code code we refer to Sec. II of 
Ref. [70]. Except for the use of sixth instead of fourth 
order spatial discretization [97 , we did not find it nec- 
essary to apply any modifications relative to the simula- 
tions presented in that work. 

The calculation of the 4-momentum in the Lean code 
is performed in accordance with the relations listed in 
Sec. |ll| The only difference is that in a BSSN code the 
four metric and its derivatives are not directly available 
but need to be expressed in terms of the 3-metric 7^^, 
the extrinsic curvature Kij as well as the gauge variables 
lapse a and shift p\ The key quantity for the calcula- 
tion of the 4-momentum is the integrand in Eq. A 
straightforward calculation gives it in terms of the canon- 
ical ADM variables 



1 



(B21) 



-P'daH^''^^, (B22) 
where K := and x •= det7~^/^ have been used for 
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convenience because they are fundamental variables in 
our BSSN implementation. 

APPENDIX C: TREATMENT OF THE EVENT 
HORIZON 

In this subsection, we summarize the numerical meth- 
ods used to find the event horizon in our pseudospectral 
simulations. 

The event horizon of the merging black holes is de- 
termined by the global structure of the spacetime, and 
thus identifying its location on any given time slice re- 
quires knowledge of the full evolution of the spacetime. 
In order to determine the event horizon surface for our 
pseudospectral evolutions, we use the "geodesic method" 
implemented by Cohen, Pfeiffer and Scheel [98], which lo- 
cates the event horizon by evolving null geodesies back- 
wards in time. This algorithm makes use of the well- 
established property that outgoing null geodesies in close 
proximity to the event horizon diverge exponentially from 
it, when followed forwards in time. Thus, these geodesies, 
when followed backwards in time, converge exponentially 
onto the event horizon, as first recognized by Libson et. 
al. [99l llQQj . One must evolve many geodesies to get a 
full picture of the horizon. In this subsection, we explain 
how these geodesies are chosen. 

The event horizon finding process can be summarized 
in three steps: 1) choosing a suitable locus of geodesies 
such that when evolved backward in time, they map out 
the event horizon, 2) evolving those geodesies backwards 
in time from a late enough time such that the spacetime 
is no longer dynamical (i.e. after the merged hole has 
rung down to its final state), and 3) determining which 
of those geodesies are on the event horizon at any given 
time. 

One property of event horizons is the formation of 
"cusps" in the course of a black hole merger. As the holes 
approach merger, generators enter the horizon through 
these cusps [101 . Therefore, as we follow the generators 
backwards in time, they will leave the horizon, and be- 
come future generators of the event horizon. This implies 
that at early times, the event horizon will be a subset of 
the surface defined by the geodesies. 

We choose our locus of geodesies from the apparent 
horizon surface at a time tend = TS.OMadm- If the co- 
ordinate location of a point on the surface is g*, we can 
expand the surface in terms of scalar spherical harmonics 
as 

L £ 

q\t,U,v) = Y, E ^Imm^mM, (CI) 

where Y^j^iu^v) are the standard spherical harmonics, 
though u and v are not the standard spherical angular 
coordinates; they are merely a conveniently chosen pa- 
rameterization of the surface. 

Individual geodesies are placed on a rectangular grid in 
(ix, v) of dimension L + 1 by 2(1/ + 1), with the u values 



chosen such that cos{u) are the roots of the Legendre 
polynomials of order L+1, and the v values are equally 
spaced in the interval [0, 27r] [98]. Therefore, there are 
N = 2(L + 1)^ geodesies on this surface; we call the value 
of L the resolution of the event horizon finding run. The 
position of the geodesic in space is given as a 3- vector 
which is evolved, along with its derivatives. The initial 
velocity of the geodesies is chosen to be the outgoing 
null normal to the surface of the apparent horizon at 

tend = 78.0MadM. 

The integration of null geodesies backwards in time is 
straightforward given the metric data from a simulation. 
Writing the position of the null geodesic as (where 

:= t), one can reexpress the geodesic equation in terms 
of coordinate time as 

= r'aprq^q' - K^H^ , (C2) 

where F^^ are the spacetime Christoffel symbols. Reex- 
pressing this as a first order system, we have 

q' = p' (C3a) 

f = r%py^p' - ri^jpV. (C3b) 

Finally, we must determine at what times some of the 
geodesies pass through the cusps of the event horizon, 
and leave the horizon (as we follow them backwards in 
time). Our most useful tool for this is the surface area el- 
ement of the event horizon. This is defined as the square 
root of the determinant of the induced metric on the 
horizon V^, where 

h = ^ dct ( lijduq\dyq{ \ 

sin^^x \7ijdvq'duq^ Jijdyq'dyq^ J ' 

^ij is the 3-metric, and the area of the event horizon is 
given by 

A{t) = JdA = J y^h{t,u,v)smudu dv. (C5) 

To determine when a null geodesic leaves the event 
horizon, (going backward in time), we note that all 
geodesies leave the event horizon surface at a cusp. In 
our pseudospectral simulations, we observe that all cusps 
on the event horizon surface are also caustics. Thus, our 
cusps may be identified by considering what happens to 
the area element of the surface at a caustic: it goes to zero 
(see §4.4 of [98] for a more thorough discussion). Thus, 
by tracking the local area element vTi of each geodesic, 
we are able to tell that it leaves the horizon at the time its 
area element approaches zero. If the local area element 
does not approach zero at any time, then it represents a 
null generator that originated on one of the two initial 
holes. We are then able to define a masking function for 
each geodesic; this function tells us if the geodesic is on 
the horizon at a given timestep. 

At this point we have located the event horizon sur- 
face at all times, and thus we may calculate the surface 
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integral of any quantity we wish on the horizon. We 
note one subtlety: the formation of cusps on the event 
horizon surface introduces a problem for taking spectral 
derivatives on the surface. Thus, for calculating deriva- 
tives of quantities on the event horizon surface, we use a 
6th order finite differencing stencil. Note that this is an 
improvement on the 2nd order stencil used in ^98J. 



APPENDIX D: 
SUPERPOSED-HARMONIC-KERR (SHK) 
INITIAL DATA 

We also present a simulation, HI in Table [T[ that is sim- 
ilar to SI except that the initial separation between the 
holes is larger and the gauge is nearly harmonic. The con- 
struction of this Superposed-Harmonic-Kerr initial data 
for this run follows that of the Superposed-Kerr-Schild 
(SI) initial data described in Appendix [Aj The differ- 
ences are as follows. 

The first difference is our choice of coordinates. In Ap- 
pendixjAj t he q uantit ies 7^ ,-, Ka, and that appear in 



Eqs. (Ao), (A6), and (All) refer to the three-metric, the 



trace of the extrinsic curvature, and the lapse function 
of the Kerr metric in Kerr-Schild coordinates. Here we 
still use Eqs. ( A5 ) , ( A6 ) , and (All), but jfj , Ka , and aa 
now refer to the three- metric, the trace of the extrinsic 
curvature, and the lapse function of the Kerr metric in 
fully harmonic coordinates, Eqs. (22)-(31), (41) and (43) 
of Ref. jlQ2j . Furthermore, the computational domain is 
excised on surfaces of constant Boyer-Lindquist radius, 
tbl, where 



1. (Dl) 



The initial coordinate separation was chosen to be 
d = 29. 73 Mad M and the Gaus sia n w eighting para m- 
eter that appears in Eqs. (|A5|), (|A6|), and ( [All ) is 



Wa = d/9. To obtain S/Ml^^ = {0, , ±0.5} we choose 
Qr = t0.261332/Madm in Eq. ( |A12| , and to make the 
total momentum vanish we chooseV^ = -0.0000582185 



Solving the XCTS equations results in initial data that 
is approximately harmonic. Harmonic coordinates satisfy 
V^Vc^^ = 0, or equivalently, := Tab^ = 0. We can 
evaluate the degree to which the harmonic gauge con- 
dition is satisfied in our initial data by examining the 
normalized magnitude of Fa'. 



,2 \ 



1/2 



\ a b / 



(D2) 



The denominator consists of the sum of squares of 
terms that must cancel to produce = 0, so that / = 1 
corresponds to complete violation of the harmonic coor- 
dinate condition. On the apparent horizons / < 0.049, 
while in the asymptotically fiat region far from the holes 
/ < 0.0083. In the regions where the Gaussians in 



Eqs. (A5), (A6) and (All) transition the XCTS free data 



from harmonic Kerr to conformally fiat we cannot ex- 
pect the data to be strongly harmonic, and we find that 
/<0.12. 

The techniques employed in the spectral evolution 
from this SHK initial data follow those used for the SKS 
initial data as described in Appendix |B 1| In particu- 
lar, the generalized harmonic gauge source function. Ha 
(Eq. 14), is constructed by demanding that Ha' remains 



in Eq. (A15). 



frozen to its value in the initial data. The evolution pro- 
ceeds in nearly harmonic gauge because of the way the 
initial data is constructed. 

Three of these HI evolutions were performed at resolu- 
tions of approximately 61^, 67^ and 72^ grid points. The 
constraints were found to be convergent. The data pre- 
sented in this paper is taken from the highest resolution 
run. 

These simulations are specifically constructed to pro- 
vide data for comparison with FN approximations, so 
we are restricted to remain in our approximately har- 
monic gauge. However, currently this gauge choice pre- 
vents us from continuing our HI evolutions beyond the 
plunge phase; we have not observed the formation of a 
common horizon. 



[1] F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005). 

[2] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlo- 

chower, Phys. Rev. Lett. 96, 111101 (2006). 
[3] J. G. Baker, J. Centrella, D.-L Choi, M. Koppitz, and 

J. van Meter, Phys. Rev. Lett. 96, 111102 (2006). 
[4] M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. 

Matthews, and H. P. Pfeiffer, Phys. Rev. D 79, 024003 

(2009). 

[5] D. Keppel, D. A. Nichols, Y. Chen, and K. S. Thorne 
(2009), submitted to Phys. Rev. D, arXiv:0902.4077. 

[6] B. Krishnan, C. O. Lousto, and Y. Zlochower, Phys. 
Rev. D76, 081501(R) (2007). 



[7] L. B. Szabados, Living. Rev. Rel. 7, 4 (2004), URL 
(citedonMay27, 2009)http : //www. livingreviews . | 
org/lrr-2004-4 
[8] C. O. Lousto and Y. Zlochower, Phys. Rev. D 79, 

064018 (2009), arXiv:0805.0159. 
[9] C. O. Lousto and Y. Zlochower, Phys. Rev. D 77, 
044028 (2008), arXiv:0708.4048. 
[10] M. Boylan-Kolchin, C.-P. Ma, and E. Quataert, Astro- 

phys. J 613, L37 (2004), astro-ph/0407488. 
[11] A. Gualandris and D. Merritt, Astrophys. J 678, 780 

(2008), arXiv:0708.0771 [astro-ph]. 
[12] S. Komossa and D. Merritt, Astrophys. J. Lett. 689, 



24 



189 (2008), arXiv:0811.1037 [astro-ph]. 
[13] P. Madau, M. J. Rees, M. Volonteri, F. Haardt, and 

S. P. Oh, Astrophys. J. 604, 484 (2004). 
[14] Z. Haiman, Astrophys. J 613, 36 (2004), astro- 

ph/0404196. 

[15] P. Madau and E. Quataert, Astrophys. J. 606, L17 
(2004). 

[16] D. Merritt, M. Milosavljevic, M. Favata, S. A. Hughes, 
and D. E. Holz, Astrophys. J. 607, L9 (2004). 

[17] M. Volonteri, Astrophys. J 663, L5 (2007), astro- 
ph/0703180. 

[18] L. Blecha and A. Loeb (2008), arXiv:0805.1420 [astro- 
ph]. 

[19] A. Loeb, Phys. Rev. Lett. 99, 041103 (2007), astro- 
ph/0703722. 

[20] S. Komossa, H. Zhou, and H. Lu, Astrophys. J. 678, 
LSI (2008), arXiv:0804.4585 [astro-ph]. 

[21] K. Menou, Z. Haiman, and B. Kocsis, New Astron. Rev. 
51, 884 (2008), arXiv:0803.3627 [astro-ph]. 

[22] M. J. Fitchett, MNRAS 203, 1049 (1983). 

[23] M. Favata, S. A. Hughes, and D. E. Holz, Astrophys. J. 
607, L5 (2004). 

[24] L. Blanchet, M. S. S. Qusailah, and C. M. WiU, Astro- 
phys. J. 635, 508 (2005), arXiv:astro-ph/0507692. 

[25] T. Damour and A. Gopakumar, Phys. Rev. D 73, 
124006 (2006). 

[26] J. G. Baker, J. Centrella, D.-L Choi, M. Koppitz, J. R. 
van Meter, and M. C. Miller, Astrophys. J. 653, L93 

(2006) , astro-ph/0603204. 

[27] J. A. Gonzalez, U. Sperhake, B. Briigmann, M. Han- 
nam, and S. Husa, Phys. Rev. Lett. 98, 091101 (2007), 
gr-qc/0702052. 

[28] F. Herrmann, L Hinder, D. Shoemaker, P. Laguna, 
and R. A. Matzner, Astrophys. J. 661, 430 (2007), gr- 
qc/0701143. 

[29] M. Koppitz, D. Pollney, C. Reisswig, L. Rezzolla, 

J. Thornburg, P. Diener, and E. Schnetter, Phys. Rev. 

Lett. 99, 041102 (2007), gr-qc/0701163. 
[30] M. Campanelli, C. O. Lousto, Y. Zlochower, and 

D. Merritt, Astrophys. J. Lett. 659, L5 (2007). 
[31] W. Tichy and P. Marronetti, Phys. Rev. D76, 

061502(R) (2007). 
[32] J. A. Gonzalez, M. D. Hannam, U. Sperhake, 

B. Briigmann, and S. Husa, Phys. Rev. Lett. 98, 231101 

(2007) , gr-qc/0702052. 

[33] M. Campanelli, C. O. Lousto, Y. Zlochower, and 
D. Merritt, Phys. Rev. Lett. 98, 231102 (2007), gr- 
qc/0702133. 

[34] J. Healy, F. Herrman, L Hinder, D. M. Schoemaker, 
P. Laguna, and R. A. Matzner, Superkicks in hyperbolic 
encounters of binary black holes (2008), 0807.3292. 

[35] L. Boyle, M. Kesden, and S. Nissanke, Phys. Rev. Lett. 
100, 151101 (2008), arXiv:0709.0299 [gr-qc]. 

[36] L. Boyle and M. Kesden, Phys. Rev. D 78, 024017 

(2008) , arXiv:0712.2819 [astro-ph]. 

[37] J. D. Schnittman and A. Buonanno (2007), astro- 
ph/0702641. 

[38] J. G. Baker, W. D. Boggs, J. Centrella, B. J. Kelly, 
S. T. McWilhams, M. C. Miller, and J. R. van Meter, 
Astrophys. J. 682, L29 (2008), arXiv:0802.0416. 

[39] W. Tichy and P. Marronetti, Phys. Rev. D 78, 
081501 (R) (2008), arXiv:0807.2985 [gr-qc]. 

[40] C. O. Lousto, M. Campanelh, and Y. Zlochower (2009), 
arXiv:0904.3541 [gr-qc]. 



[41] J. G. Baker, W. D. Boggs, J. Centrella, B. J. Kelly, 
S. T. McWilliams, M. C. Miller, and J. R. van Meter, 
Astrophys. J. 668, 1140 (2007), astro-ph/0702390. 

[42] J. A. Gonzalez, U. Sperhake, and B. Briigmann, Phys. 
Rev. D 79, 124006 (2009), arXiv:0811.3952 [gr-qc]. 

[43] L. Rezzolla, Class. Quant. Grav. 26, 094023 (2009), 
arXiv:0812.2325 [gr-qc]. 

[44] J. D. Schnittman, A. Buonanno, J. R. van Meter, J. G. 
Baker, W. D. Boggs, J. Centrella, B. J. Kelly, and S. T. 
McWilliams, Phys. Rev. D 77, 044031 (2008). 

[45] S. H. Miller and R. Matzner, Gen. Rel. Grav. 41, 525 
(2009), arXiv:0807.3028. 

[46] E. Racine, A. Buonanno, and L. E. Kidder (2008), ac- 
cepted for publication in Phys. Rev. D, arXiv:0812.4413. 

[47] Y. Mino and J. Brink, Phys. Rev. D 78, 124015 (2008), 
arXiv:0809.2814. 

[48] F. Pretorius (2007), arXiv:0710.1338 [gr-qc]. 

[49] L. D. Landau and E. M. Lifshitz, Classical Theory of 
Fields (Addison Wesley, Reading, MA, 1962), 2nd ed. 

[50] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Grav- 
itation (Freeman, New York, New York, 1973). 

[51] S. V. Babak and L. P. Grishchuk, Phys. Rev. D 61, 
024038 (1999). 

[52] G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, Phys. 

Rev. D 78, 084017 (2008). 
[53] R. M. Wald, General Relativity (University of Chicago 

Press, 1984). 

[54] R. Arnowitt, S. Deser, and C. W. Misner, in Gravita- 
tion: An Introduction to Current Research, edited by 
L. Witten (Wiley, New York, 1962), pp. 227-265, gr- 
qc/0405109. 

[55] J. W. York, Jr., in Sources of Gravitational Radiation, 

edited by L. L. Smarr (Cambridge University Press, 

Cambridge, England, 1979), pp. 83-126. 
[56] H. Friedrich, Commun. Math. Phys. 100, 525 (1985). 
[57] D. Garfinkle, Phys. Rev. D 65, 044029 (2002). 
[58] F. Pretorius, Class. Quantum Grav. 22, 425 (2005). 
[59] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, and 

O. Rinne, Class. Quantum Grav. 23, S447 (2006). 
[60] M. Ruiz, R. Takahashi, M. Alcubierre, and D. Nunez, 

Gen. Relativ. Gravit. 40, 1705 (2008), arXiv:0707.4654. 
[61] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. 

D 56, 3405 (1997). 
[62] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, 

D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 

67, 084023 (2003). 
[63] S. Brandt and B. Briigmann, Phys. Rev. Lett. 78, 3606 

(1997). 

[64] M. Ansorg, B. Brii gmann, and W. Tichy, Phys. Rev. D 

70, 064011 (2004), gr-qc/0404056. 
[65] J. M. Bowen and J. W. York, Jr., Phys. Rev. D 21, 2047 

(1980). 

[66] M. J. Berger and J. Ohger, J. Comput. Phys. 53, 484 
(1984). 

[67] E. Schnetter, S. H. Hawley, and L Hawke, Class. Quan- 
tum Grav. 21, 1465 (2004). 

[68] Carpet - Adaptive Mesh Refinement for the Cactus 
Framework, http:/ /www. carpetcode.org. 

[69] M. Alcubierre, B. Briigmann, P. Diener, M. Koppitz, 
D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 
67, 084023 (2003). 

[70] U. Sperhake, Phys. Rev. D 76, 104015 (2007), gr- 
qc/0606079. 

[71] G. Faye, L. Blanchet, and A. Buonanno, Phys. Rev. D 



74, 104033 (2006). 
[72] H. Tagoshi, A. Ohashi, and B. J. Owen, Phys. Rev. D 

63, 044006 (2001). 
[73] G. B. Cook, Phys. Rev. D 65, 084003 (2002). 
[74] G. B. Cook and H. P. Pfeiffer, Phys. Rev. D 70, 104016 

(2004). 

[75] M. Caudill, G. B. Cook, J. D. Grigsby, and H. P. Pfeif- 
fer, Phys. Rev. D 74, 064011 (2006). 

[76] E. Gourgoulhon, P. Grandclement, and S. Bonazzola, 
Phys. Rev. D 65, 044020 (2002). 

[77] P. Grandclement, E. Gourgoulhon, and S. Bonazzola, 
Phys. Rev. D 65, 044021 (2002). 

[78] H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. 
Teukolsky, Comput. Phys. Commun. 152, 253 (2003). 

[79] M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder, 
O. Rinne, and S. A. Teukolsky, Phys. Rev. D 74, 104006 
(2006). 

[80] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. 
Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, 
Phys. Rev. D 76, 124038 (2007). 

[81] O. Rinne, Class. Quantum Grav. 23, 6275 (2006). 

[82] O. Rinne, L. Lindblom, and M. A. Scheel, Class. Quan- 
tum Grav. 24, 4053 (2007). 

[83] J. M. Stewart, Class. Quantum Grav. 15, 2865 (1998). 

[84] H. Friedrich and G. Nagy, Commun. Math. Phys. 201, 
619 (1999). 

[85] J. M. Bardeen and L. T. Buchman, Phys. Rev. D 65, 
064037 (2002). 

[86] B. Szilagyi, B. Schmidt, and J. Winicour, Phys. Rev. D 

65, 064015 (2002). 
[87] G. Calabrese, J. Pullin, O. Reula, O. Sarbach, and 

M. Tiglio, cmp 240, 377 (2003), gr-qc/0209017. 
[88] B. Szilagyi and J. Winicour, Phys. Rev. D 68, 

041501(R) (2003). 
[89] L. E. Kidder, L. Lindblom, M. A. Scheel, L. T. Buch- 
man, and H. P. Pfeiffer, Phys. Rev. D 71, 064020 (2005). 
[90] L. T. Buchman and O. C. A. Sarbach, Class. Quantum 

Grav. 23, 6709 (2006). 
[91] L. T. Buchman and O. C. A. Sarbach, Class. Quantum 

Grav. 24, S307 (2007). 
[92] U. Sperhake, B. Briigmann, J. Gonzalez, M. Hannam, 

and S. Husa, in Proceedings of the eleventh Marcel 

Grossmann Meeting (2007), 0705.2035vl. 
[93] The cactus computational toolkit^ 

http:/ /www. cactuscode.org. 
[94] J. Thornburg, Phys. Rev. D 54, 4899 (1996), gr- 

qc/9508014. 

[95] J. Thornburg, Class. Quantum Grav. 21, 743 (2004), 
gr-qc/0306056. 

[96] J. R. van Meter, J. G. Baker, M. Koppitz, and D.-L 
Choi, Phys. Rev. D 73, 124011 (2006), gr-qc/0605030. 

[97] S. Husa, J. A. Gonzalez, M. Hannam, B. Briigmann, and 
U. Sperhake, Class. Quantum Grav. 25, 105006 (2008). 

[98] M. A. S. Michael Cohen, Harald P. Pfeiffer, Class. Quan- 
tum Grav. 26, 035005 (2009), arXiv:0809.2628. 

[99] J. Libson, J. Masso, , E. Seidel, W.-M. Suen, and 
P. Walker, Phys. Rev. D53, 4335 (1996). 
[100] P. . Anninos, D. . Bernstein, S. . Brandt, J. . Libson, 
J. . Masso, E. . Seidel, L. . Smarr, W.-M. . Suen, and 
P. . Walker, Phys. Rev. Lett. 74, 630 (1995). 
[101] S. L. Shapiro, S. A. Teukolsky, and J. Winicour, Phys. 

Rev. D 52, 6982 (1995). 
[102] G. B. Cook and M. A. Scheel, Phys. Rev. D 56, 4775 
(1997). 



